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Pre  f  ace 


The  purpose  of  this  study  was  to  determine  how  mass 
motion  on  a  gravity  gradient  stabilized  space  station  would 
alter  the  classical  Lagrange  stability  region  described  by 
(Debra  and  Delp,  1961).  The  addition  of  mass  motion  to  an 
orbiting  space  station  introduces  time  varying  moments  of 
inertia  into  the  equations  of  attitude  motion.  When  the 
time  variance  is  periodic,  it  is  a  source  of  parametric 
excitation. 

The  phenomenon  of  a  parametrically  excited  system  is 
also  seen  when  eccentricity  is  introduced  into  the  orbit  of 
a  gravity-gradient  stabilized  rigid  body  as  shown  by 
( Breakwel 1  and  Pringle,  1965).  This  study  obtained  similar 
results  and  shows  that  if  masses  move  across  the  surface  of 
a  future  space  station,  the  gravity  gradient  torque  will 
cause  unstable  attitude  motion  for  certain  configurations  in 
the  classical  Lagrange  region. 

I  would  like  to  thank  my  thesis  advisor,  Dr.  Curtis  H. 
Spenny  of  the  Department  of  Aeronautics  and  Astronautics, 
for  his  valued  assistance  in  guiding  me  through  this  study. 

I  would  also  like  to  thank  Captain  James  Planeaux,  also  of 
the  Department  of  Aeronautics  and  Astronautics,  for  his 
critical  review  of  the  final  product.  Finally,  I  would  like 
to  thank  my  wife  Beverly  for  her  patience  and  support; 
without  it,  this  study  would  not  have  been  possible. 
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Abst  rac t 


The  future  space  station  will  contain  a  mechanism  that 
transports  mass  across  large  distances  of  its  surface. 
Accordingly,  this  study  will  derive  the  equations  of 
attitude  motion  for  a  gravity-gradient  stabilized  space 
station  whose  moments  of  inertia  are  varying  with  time.  The 
equations  are  linearized,  after  which  Hill’s  Equation  is 
used  to  determine  pitch  stability,  while  the  Method  of 
Multiple  Scales  is  used  to  determine  the  roll/yaw  stability 
of  the  system.  Results  show  that  for  certain  frequencies  of 
mass  motion,  attitude  motion  can  grow  unboundedly  with  time. 
Consequently,  the  shape  of  the  classical  Lagrange  stability 


region  is  altered. 


AN  ANALYSIS  OF  THE  ATTITUDE  STABILITY 
OF  A  SPACE  STATION  SUBJECT  TO 
PARAMETRIC  EXCITATION  OF 
PERIODIC  MASS  MOTION 


I .  Introduction 

The  proposed  U.S.  space  station  will  be  larger  than  any 
spacecraft  previously  built  by  man.  As  the  station  grows,  it 
will  become  necessary  to  move  maintenance  and  modification 
equipment  over  large  distances  of  the  space  station  surface. 
Hence,  there  is  a  need  for  a  trolley  apparatus  that  can 
transport  large  masses  from  section  to  section.  This  mass 
movement  will  affect  the  attitude  stability  of  the  space 
station.  Research,  therefore,  is  needed  to  determine  the 
magnitude  of  this  destabilizing  motion  and  how  it  can  be 
minimized,  in  order  to  design  an  adequate  control  system. 

Objective 

The  objective  of  this  study  is  to  examine  what  happens 
to  the  pitch,  roll,  and  yaw  stability  of  a  gravity-gradient 
stabilized  space  station  when  a  trolley  apparatus  operates 
along  a  defined  path.  This  study  expands  on  the  work  John 
Chan  (Chan,  1986)  did  with  the  pitch  equation.  The  analysis 
continues  by  investigating  the  behavior  of  the  coupled 
roll/yaw  equations  of  attitude  motion  in  order  to  define 
comprehensive  conditions  where  mass  movement  will  not  lead 
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to  large  oscillations  or  to  totally  unstable  attitude 


mot  ion . 

Methodology 

The  equations  of  motion  are  derived  for  a  gravity- 
gradient  stabilized  space  station  whose  moments  of  inertia 
are  varying  with  time.  The  equations  are  then  linearized, 
after  which  Hill’s  Equation  and  the  Method  of  Multiple 
Scales  are  used  to  determine  the  pitch,  roll,  and  yaw- 
stability  of  the  system.  Stability  boundaries  are 
determined  so  that  the  effect  of  mass  motion  on  the 
classical  Lagrange  gravity-gradient  stability  region  is 
observed.  The  results  are  compared  to  the  solutions 
obtained  from  numerical  integration  of  the  nonlinear 
equations  for  verification. 

Scope 

This  study  is  restricted  by  the  following  conditions: 
(1)  The  only  external  torques  acting  on  the  spacecraft  are 
gravitational  environmental  torques;  (2)  The  spacecraft  is 
small  enough  so  that  attitude  motions  have  no  significant 
effect  on  the  orbital  motion  of  the  center  of  mass;  (3)  The 
orbit  is  circular. 


I I .  Background  Theory 

The  purpose  of  this  Chapter  is  to  provide  some 
background  theory  on  gravity-gradient  stabilization  and 
torques  on  spacecraft  due  to  mass  motion.  An  overview  is 
provided  on  the  analytical  methods  used  to  determine  the 
stability  of  the  linearized  equations  of  motion.  Specifics 
on  each  of  the  methods,  however,  are  provided  in  Chapters  V 
and  V I  . 

Spacecraf t  Torques 

The  purpose  of  this  section  is  to  survey  specific 
spacecraft  torques  and  describe  their  origins  and 
characteristics.  (Hughes,  1986:  232)  points  out  that  the 
most  striking  characteristic  of  spacecraft  torques  is  their 
minuteness.  In  terras  of  familiar  terrestrial  experience, 
they  are  intuitively  negligible.  On  closer  examination 
however,  there  are  no  "large"  torques  in  space,  and  hence 
minor  influences  play  major  roles  in  governing  the  attitude 
dynamics  of  spacecraft.  To  help  distinguish  between  the 
many  different  torques  acting  on  spacecraft,  they  are 
divided  into  two  categories:  internal  and  external. 

Internal  Torques 

Internal  torques  are  usually  defined  as  torques 


exerted  on  the  main  body  of  a  spacecraft  by  internal  moving 


parts.  Internal  moving  parts  include  reaction  wheels, 
flexible  booms  or  solar  arrays,  scanning  or  rastering 
instruments,  liquids  inside  partially  filled  tanks,  and 
moving  crew  members.  For  the  purposes  of  this  study, 
internal  torques  will  also  include  torques  caused  by  the 
motion  of  mass  along  the  surface  of  a  spacecraft.  This  is 
an  obvious  extension  because  all  internal  torques  on 
spacecraft  are  in  a  sense  se 1 f -generated . 

The  effect  of  mass  motion  on  the  stability  of 
spacecraft  has  been  studied  before.  Research  and  experience 
have  shown  that  disturbances  due  to  mass  motion  are  a 
function  of  the  amplitude  of  the  motion  and  the  ratio  of  the 
moving  mass  to  the  spacecraft’s  moment  of  inertia. 

(Thomson,  1985:  1087)  showed  analytically  that  if  crew 

members  on  a  space  station  moved  periodically,  they  could 
rock  a  station  and  cause  it  to  tumble.  In  addition,  (Chubb, 
1975:  277)  used  flight  data  from  Skylab  to  show  that  while 
crew  members  were  jogging,  the  station  experienced 
substantial  angular  deviation  from  its  nominal  attitude. 

It  should  be  pointed  out,  however,  that  the  crew  motion 
discussed  in  the  previous  articles  was  over  small  distances 
and  that  the  ratios  of  the  crew  members’  inertias  to  the 
space  stations’  moment  of  inertias  were  small.  Research, 
therefore,  is  required  to  examine  the  movement  of  relatively 
large  masses  over  large  distances. 
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External  Torques 


External  torques  arise  through  the  interaction  of  a 
vehicle  with  its  environment.  They  include  gravitational, 
aerodynamic,  radiation,  and  magnetic  torques.  This  study 
will  concentrate  specifically  on  the  gravitational 
environmental  torque. 

The  Earth’s  g rav i tat ional  field  is  one  of  the  dominant 
sources  of  attitude  disturbance  torques  on  a  spacecraft  in 
low-earth  orbit.  Any  nonsymmetrical  object  of  finite 
dimensions  in  orbit  is  subject  to  a  gravitational  torque 
because  of  the  variation  in  the  Earth’s  gravitational  force 
over  the  object.  This  variation  is  caused  by  the  inverse 
square  gravitational  force  field  that  surrounds  the  earth. 
The  result  is  a  gravity-gradient  torque  which  will  try  to 
align  a  spacecraft’s  minimum  moment  of  inertia  axis  with  the 
local  vertical  (radius  vector  drawn  from  center  of  the 
earth ) ( Wertz,  1985:  566). 

The  principle  of  gravity-gradient  stabilization  can  be 
explained  simply  by  considering  the  attitude  motion  of  the 
dumbbell  satellite  in  Fig.  1,  where  the  orbital  frame  of 
reference  is  o^  o2,  and  o3  with  the  origin  at  the  center  of 
mass  of  the  body,  the  axis  ot  directed  away  from  the  center 
of  the  earth,  axis  o2  tangent  to  the  orbit  in  the  direction 


of  motion,  and  axis  o  normal  to  the  orbit  plane.  The  body 

A  A  A 

principal  axes  are  represented  by  b1 ,  b2 ,  and  b3 .  A 
deflection  away  from  the  local  vertical  causes  a  restoring 
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torque  to  be  generated  by  the  imbalance  of  forces  acting  on 
the  spherical  masses.  The  centrifugal  force  acting  on  the 
sphere  furthest  from  the  earth  is  greater  than  the 
gravitational  force  on  it  because  these  two  forces  are  equal 
only  at  the  center  of  mass.  Since  the  oposite  is  true  of 
the  mass  closet  to  the  earth,  a  net  torque  is  created  which 
forces  the  masses  toward  a  local  vertical  orientation 
(Kaplan,  1976:  199). 

Gravity-gradient  stabilization  works  not  only  for 
dumbbell  satellites  but  for  any  spacecraft  whose 
configuration  lies  within  the  regions  shown  in  Fig.  2.  For  a 
rigid  body  in  a  circular  orbit,  the  Lagrange  region  is 
statically  stable  whereas  the  Delp  region  is  gyrically 
stable.  Any  inertia  shape  that  lies  outside  of  these 
regions  would  result  in  unstable  motion.  Moreover,  when 
flexibility  effects  are  considered  with  a  quasi-rigid  body, 
the  Delp  region  loses  its  gyrically  stable  character  while 
the  Lagrange  region  becomes  asymptotically  three-axis  stable 
(Hughes,  1986:  318).  At  the  same  time,  Chan  (20)  points  out 

A 

that  having  mass  motion  along  the  body  axis  b^  (see  Fig.  1) 
of  a  space  station  whose  inertia  shape  is  in  the  Delp  region 
would  correspond  to  mass  motion  along  the  shortest  axis  of 
the  station  inertia  ellipsoid.  Therefore,  the  Delp  region 
will  not  be  considered  here. 
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Fig.  2. 

Attitude  Stability  and  Natural  Frequencies  of  a  Rigid  Body 
in  a  Circular  Orbit  Plotted  Versus  ri  and  r* 
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Anal vt ical  Methods 


A 

When  an  elevator  that  moves  along  the  b  axis  is  added 
to  a  dumbbell  satellite  (Fig.  3),  the  satellite  will  "bob" 
up  and  down  along  this  axis  in  order  to  maintain  the  center 
of  mass  orbital  path.  If  the  satellite  is  deflected  from 
the  local  vertical,  however,  periodic  mass  movement  will 
exaggerate  the  torque  caused  by  the  gravity-gradient  and 
result  in  large  oscillations  in  the  satellite’s  attitude. 
Tins  phenomenon  is  known  as  parametric  excitation  and 
results  when  the  coefficients  of  the  equations  of  attitude 
motion  vary  periodically  with  time. 

In  an  attempt  to  find  an  analytical  solution  for  these 
equations,  they  must  first  be  linearized.  The  results, 
detailed  in  Chapter  IV,  show  linearized  pitch,  roll,  and  yaw 
equations  where  the  roll  and  yaw  equations  are  coupled. 

Many  theories  exist  that  determine  the  behavior  of  systems 
governed  by  linear  ordinary  differential  equations  with 
periodic  coefficients.  This  study,  however,  will  focus  on 
Hill’s  equation  and  a  generalized  form  of  Hill’s  equation 
called  Mathieu’s  equation  that  (Chan,  1986)  used  to  conduct 
the  pitch  analysis.  In  addition,  the  Method  of  Multiple 
Scales  is  used  to  analyze  the  coupled  roll/yaw  equations. 
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III.  Per i vat  ion  o f  the  Eguat ions  o f  Motion 


The  equations  of  attitude  motion  of  a  gravity-gradient 
stabilized  spacecraft  of  arbitrary  shape,  for  the  conditions 
specified  in  Chapter  I,  are  well  known  and  can  be  seen  in 
(Kaplan,  1976:  201).  A  difference  occurs  in  this  study, 
however,  when  a  moving  mass  is  added.  This  difference 
results  in  the  principal  moments  of  inertia  of  the 
spacecraft  becoming  a  function  of  time. 

The  time  varying  moments  of  inertia  were  taken  into 


account  b \ 

(  Chan  , 

1986  ) 

and  resulted  in  three  Euler 

moment 

equat ions 

V,  + 

V,  + 

(I3- 

Ws  = 

( 3M/R5 ) ( I3- 

(  1  ) 

+ 

I2W2  + 

(Ir 

Vv'i  = 

(  3M/R5  )(I,- 

I3)R,R3 

(  2  ) 

v3  + 

V3  + 

(I2~ 

I1,WlW2  = 

(3M/R5)(I2- 

vb,r2 

(  3  ) 

Eqs .  (1),  (2),  and  (3)  differ  from  standard  Euler  equations 

for  a  gravity-gradient  stabilized  body  (Kaplan,  1976:  201) 
because  of  the  Iw  terms  which  permit  the  introduction  of 
mass  movement  along  an  arbitrary  axis.  Their  derivation  is 
summarized  below  for  completeness. 

The  formulation  begins  with  an  inverse  square 
formulation  of  gravitational  torques  given  in  vector  form  by 
Wertz  ( 567  ) 

M  =  (3JJ/R5)  [R  *  (I-BH 


(4) 


where  R  is  the  radius  vector  from  earth  center  to  the  mass 
center  of  the  spacecraft,  and  /J=GM  is  the  earth’s 
gravitational  constant.  Since  it  is  stated  in  Chapter  I 


that  there 

is  no  attitude/orbit 

coupling,  Eq.  (4)  in 

body 

component 

form 

becomes 

M1  =  ( 3M/R5  )[  ( I3- 

I2)R2R3] 

(  5c  ) 

M2  =  (3M/R6)[ ( Ij- 

(  5b) 

M3  =  ( 3M/R5 )( < I2- 

Ii)RiR21 

(  5c  ) 

Recall 

that 

the  reaction  of 

this  gravity-gradient 

torque 

(Mi  about  a  rigid  body’s  center  of  mass  is  equal  to  the  time 
rate  of  change  of  the  body’s  angular  momentum  (h)  with 
respect  to  an  inertial  frame 

M  ■  3t(h>  •  161 

Eq .  ( 6 )  becomes 

M  =  Tj-tl  ’w8]  +  IwB  x  [I-'wB]  (7) 

at  — 

since 

h  =  I  ■  1  wB  (  8  ) 

I  B 

where  I  is  the  inertia  dyadic  of  the  space  station  and  w 
is  the  angular  velocity  of  the  body  frame  with  respect  to 


the  inertial  frame. 


Because  the  origin  of  the  body  reference  frame  b  is 


located  at  the  geometric  center  of  the  space  station, 
motion  will  cause  it  to  translate  in  space  relative  to 
center  of  mass.  Therefore,  the  first  term  in  Eq.  (7) 
becomes 


Bd  r  T  1  B  .  T  Bd  r  1  B  t  I  B  ®d  r  T  i 

w  ]  =  !•—[  w  J  +  w  •  dtHl 


or  in  matrix  form, 


I  B 
w 


'  *1 

0 

0 

'  .  ' 

w 

1 

' 

0 

0 

C  ' 

w 

1 

0 

J2 

0 

< 

w„ 

>  + 

0 

*2 

0 

< 

u2 

0 

0 

1  3- 

1  3  J 

0 

0 

J3- 

.  w3. 

11  *11 
:2W2  +  i2K2  [ 
.  I3K3  *  i3W3  J 


In  addition,  the  second  term  of  Eq.  (7)  becomes 


or 


ib  .  _  I  B  - 
w  x  [  I  '  w  ]  = 


I  B  .IB, 
W  x  [  I  -  W  1 


A 

A 

A 

b1 

b  2 

b  3 

W1 

W2 

W3 

I1W1 

I2W 

2  J3W 

r 

<V 

1  2 

)W2W3 

-  (I,- 

J3 

)W1W3 

l  <V 

*1 

)W1W2 

mass 

the 


(  9a  ) 


(  9b) 


(  9c  ) 


(  10a  ) 


(10b) 


Substituting  Eqs.  (5),  (9c),  and  (10b)  into  Eq .  (7)  will 
yield  Eqs.  (1),  (2),  and  (3). 


As  pointed  out  by  (Chan,  1986:  6)  the  orientation  of  the 

body  cannot  be  determined  from  Eqs.  (1),  (2),  and  (3). 

Therefore,  relationships  are  needed  between  these  Euler 
moment  equations  and  angles  specifying  the  orientation  of 
the  body.  One  way  of  doing  this  is  with  Euler  angles. 

The  direction  cosine  matrix  using  a  3-1-2  sequence  of 

A 

A 

rotations  between  the  o  frame  and  the  b  frame  is 


A  3 1  z  <  p  ’  y  ’  r  >  =  A2lr)  A1(y)  A3(p) 


c  rep- sy  s  rsp 
-cysp 

srep+sversp 


ersp+sysrep  -cysr 
eyep  sy 

srsp-sycrcp  cycr 


(11) 


where  the  Euler  angles  p,  y,  and  r  represent 

P  *  pitch  angle  »  body  rotation  about  c>3 
y  *  yaw  angle  *  body  rotation  about 
r  *  roll  angle  ■  body  rotation  about  o2 
and  c  =  cosine  and  s  =  sine. 


The  expressions  for  the  rotation  angles  in  terms  of  the 
e 1 emen t s  of  the  direction  cosine  matrix  are 


-arctan 

(A2,/A2z' 

(12a) 

arc  sin 

<A23> 

(12b) 

-arctan 

(A,3/A33' 

(12c) 

The  angular  velocity  of  the  body  frame  with  respect  to 


the  orbit  frame  is 


-  p  coslyl  sin(r)  +  y  coslr) 
P  s  i  n  (  y  )  +  r 

p  cos(y)  coslr)  +  v  sin(r) 


f  0  B  i 

w 

1 

0  B 

W 

*  -  < 

0  B 

W  „ 

V  3  J 

V 

(13) 


Since  we  need  the  angular  velocity  of  the  body  frame  with 


i  B 


respect  to  the  inertial  frame  w 


IB  10  OB 
w  =  w  +  w 
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where  'w°  =  No^  (N  *  scalar  eccentric  orbital  rate). 
Therefore,  Eq .  (14)  becomes 


K1 

►  - 

A11  A 1 2  A 1 3 
A21  A22  A23 

< 

r 

O  O 

>  +  - 

-pcos ( y ) s in ( r )  +  ycos(r) 

psin( y  )  +  r 

W 

V  3j 

-A31  A32  A  33" 

l  N  J 

pcos(y)coslr)  +  ysin(r) 

or  , 


w  =  A 1  ,(N  -  pc  os  (  y  )  sin(  r  )  +  vcos(r)  (16a) 

W2  =  ‘^23^  +  psin(y)  +  r  (16b) 

W3  =  a33N  +  pcos ( y ) cos ( r  )  +  ysin(r)  (16c) 


The  Euler  angle  rates  can  now  be  solved  in  terms  of  the 
angular  velocity  and  Euler  angles: 


[(w3-  A33N 

)  cos ( r ) 

-  <w  -  A13N)  sin(r)]  sec(y) 

(17a 

(wt-  A13N) 

cos  (  r  )  + 

(w3-  A33N)  sin(r) 

(  17b 

(<w3-  A33N) 

cos ( r ) 

-  (w  -  AN)  sin(r)]  tan(y) 

1  i  3 

(17c 

+  (w2-  A23N) 

A  computer  simulation  of  a  gravity-gradient  stabilized 


space  station  with  mass  movement  can  now  be  generated  by 
solving  Eqs .  (1),  (2),  and  (3)  simultaneously  with  Eqs . 

(17a),  (17b),  (17c)  for  given  values  of  i. 

Although  the  geometric  significance  of  Euler  angles  is 
more  apparent,  computationally  they  are  not  as  efficient  as 
Euler’s  symmetric  parameters  q^  q2 ,  q3 ,  and  q^  .  Appendix  A 
details  how  Euler  symmetric  parameters  (known  as 
quaternions)  are  used  instead  of  Eqs.  (17)  in  the  nonlinear 
program  of  this  study. 
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IV .  Linearization 


The  Iw  terms  in  the  Euler  moment  Eqs .  (1),  (2),  and  (3) 

of  Chapter  III  permit  the  introduction  of  mass  movement 

along  an  arbitrary  axis.  If  a  trolley  apparatus  is 

A 

constrained  to  move  along  the  body  axis  ( b 1 ) ,  then  the  term 
I  vanishes.  This  leaves  us  with 

V’l  +  <V  VW2W3  =  (3^/R5)  (I3"  I2>R2R3  (18a) 


V2 

+  i2w2  + 

<V  V'S'S 

=  (3m/R5) 

<V 

1  3 

(  18b) 

X3W3 

+  i3W3  + 

d2-  I,)w1w2 

=  (3M/R5) 

(I2- 

ir,r2 

(  18c  ) 

For 

small  p, 

y,  and  r,  the 

direction 

cosine 

matrix 

becomes 


A  (p,y,r) 


1  p  -r 
-P  1  y 
r  -y  1 


(  19  ) 


In  order  to  get  the  components  of  the  radius  vector  R,  we 
can  write 


r  > 

R, 

l 

R 

>  ~ 

2 

l  r3J 

L 

P 

1 

-y 


(20) 


or  , 


In  addition, 


f  \ 

R, 

f 

R 

R2 

►  =  . 

-Rp 

l  r3J 

.  Rr  , 

the  angular  velocity  of 


(21  ) 

the  orbit  frame  with 
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Di f  f erent 

iating 

Eqs  . 

(  24  ) 

with 

respect  to 

time 

in  the  body 

frame  to 

obtain 

*1 

=  y 

-  (  rN 

+  rN) 

(  25a  ) 

W2 

=  r 

+  (  yN 

+  yN ) 

(  25b) 

W3 

=  P 

+  N 

(25c) 

For  small 

orbit 

eccentric 

ity  e , 

let  N  and 

R  be 

approximated 

by 

N  = 

n  [  1 

+  2ecos ( nt  )  ] 

p 

(26) 

R  = 

a  [  1 

-  ecos ( nt  ) ] 
p 

(27  ) 

where  n  = 

(M/a3) 

.  5  . 

IS 

the 

mean  motion,  and 

t  is  the  time 

n 

p 


since  periapsis  passage.  Differentiating  Eq.  (26)  with 
respect  to  time  leads  to 


18 


I 


N  =  -2n2esin(nt  )  (28) 

p 

After  substituting  Eqs.  (215,  (24),  and  (25)  into  Eqs.  (18), 

the  Euler  Moment  equations  become 

I  j  (  y-rN-Nr )+(  1 3~1 2  )  (  r+y.N  )  (  p  +  N  )  =  ( 3p/R5  )  (  1 3~ 1 2  )  (  -pR)  (  rR  )  (  29  ) 

I2<  r  +  yN  +  Ny )  +  i 2 ( r+yN )  +  ( 1 1~1 3 ) ( p+N ) ( y-rN ) 

=  (3^/R5)(I1-I3)(rR)(R)  (30) 

I 3 ( p+N )  +  i 3 ( P  +  N )  +  ( 1 2~1 1 ) ( y-rN ) ( r+yN  ) 

=  (3m/R5)(I2-I1)(R)(-PR)  (31) 

Substitute  N  and  N  from  Eqs  (26)  and  (28)  into  the  above 
equations  and  delete  the  products  of  quantities  assumed  to 

.  3 

be  small,  such  as  py  and  pe .  Finally,  substitute  H/R  for 
N2  and  Eqs.  (29),  (30),  and  (31)  become 

1  iy  +  ( 1 3~ 1 2~1 1 )nr  +  (I3-I2)n2y  =  0  (32) 

1 2 r  +  (I2  +  I1-I3)ny  -  4(I1-I;j)n2r  +  i2(r+yn)  =  0  (33) 

I _ [ p-2n2esin( nt  )]  +  3n2(I  -I,)p  +  i_(p+n)  =  0  (34) 

It  should  be  noted  that  the  pitch  equation,  Eq.  (34),  is 
uncoupled  from  the  coupled  yaw  and  roll  equations,  Eqs.  (32) 
and  (33).  In  addition,  as  stated  in  (Chan,  1986:  13),  these 
equations  are  nearly  identical  to  those  derived  by  Kaplan 
(204)  for  a  gravity-gradient  stabilized  rigid  body.  The 
only  difference  is  the  existence  of  the  I  terms  due  to  mass 


motion . 


The  objective  of  this  chapter  is  to  obtain  an  analytical 
solution  to  the  linearized  pitch  equation  using  Hill’s 
Method.  The  results  are  compared  to  those  obtained  using 
Mathieu’s  equation  (Chan,  1986)  and  to  the  numerical 
solutions  of  the  nonlinear  system. 

The  previous  work  done  by  Chan  on  this  subject  used 
Mathieu’s  equation  to  determine  the  inertia  shapes 
in  the  Lagrange  stability  region  that  would  not  remain  stable 
when  the  system  was  parametrically  excited  by  mass  movement. 
This  chapter  begins  by  reviewing  the  theory  and  derivations 
that  lead  to  the  use  of  Mathieu’s  and  Hill’s  equations  and 
then  recaps  the  results  Chan  realized.  The  remaining  work 
concentrates  on  using  Hill’s  equation  to  find  unstable 
configurations  and  comparing  the  results  to  those  of  (Chan, 
1986  )  . 


Theory /Deri vat  ions 


The  pitch  equation  was  shown  to  be  decoupled  from  the 
roll  and  yaw  equations  in  the  previous  chapter.  After 
rearranging  Eq .  (34) 


2n  esin(  nt  ) 

p 


(35) 


we  should  note  that  what  is  left  is  a  linear  ordinary 
differential  equation  of  the  form 


(  36  ) 


p  +  2P(t)p  +  R2(t)p  =  Q(t) 

where  the  coefficients  are  time  varying  and  defined  as 


follows : 

P( t )  =  .5  ( i3/  I3  )  ( 37  ) 

R2 ( t )  =  3n2[(I2-I1 )/  I3]  (38) 

Q(t)  =  — ( i _ /  I_)n  +  2n2esin(t  )  (39) 

3  3  p 

By  introducing  the  variable  transformation 

p  =  z  e  ( 40  ) 

Chan  (15)  showed  that  Eq .  (36)  becomes 

z  +  [R2-  P2-  P]z  =  Qe^Pdt  (41) 


In  addition,  if  the  mass  motion  is  chosen  to  be  sinusoidal, 
then  the  coefficient  of  z  is  a  periodic  function  of  time  and 
can  be  expanded  to 

[R2-P-P]  =  [aT+  1 6q  ^  cos  2T  +  16q2cos4T 

+  16q3cos6T  +  16q^cos8T]  (42) 

where  Appendix  B  details  this  expansion  as  well  as 
demonstrates  that  the  final  form  of  the  transformed, 
homogeneous  pitch  equation  is  that  of  Hill’s  equation 

z  +  [6  +  2)  e,cos(2iT)]z  =  0  (43) 

0  1 

where  0  =  a_  and  0  =  8q.. 

0  T  i  i 

The  general  solution  to  Hill’s  equation  is 

z  =  cteMT$  (T,a)  +  c2e'MT<f>  (  T , -a )  (44) 

where  the  stability  of  the  system  will  be  determined  by  the 


characteristic  exponent  p  (Hayashi,  1964:  341).  Meirovitch 
(282)  shows  that  the  determination  of  (J  is  not  an  easy 
matter.  However,  since  p  depends  on  the  parameters  0Qand 
0  ,  it  is  possible  to  represent  solutions  to  Eq.  (44)  by  a 
parameter  space  defined  by  these  quantities.  Meirovitch 
(282)  states  that  "the  parameter  space  can  be  divided  into 
regions  of  bounded  and  unbounded  motions  with  periodic- 
motions  providing  the  surfaces  separating  these  regions." 

If  the  parameter  space,  however,  is  truncated  to  the  two 
parameters  and  ,  the  boundary  surfaces  become  boundary 
curves.  This  lead  Chan  (18)  to  use  a  specialized  form  of 
Hill’s  equation  called  Mathieu’s  equation 

z  +  (aT+  16qicos2T)z  =  0  (45 

to  determine  the  stability  of  the  pitch  equation.  The 
result  is  Fig.  4,  which  shows  regions  of  bounded  and 
unbounded  motion  as  a  function  of  the  parameters  aT  and  q^. 
This  type  of  representation  is  referred  to  as  a  Strutt 
diagram . 


Summary  of  Previous  Resul ts 


In  order  to  simplify  the  analysis,  Chan  needed  to 
determine  which  regions  of  instability  are  of  concern. 
Appendix  B  shows  that  the  parameters  aT  and  q1  are  a 


function  of  the  variables  r  ,  r  ,  r  and  (n/w)  where  they 

2  I1  I 

,  ,  ML  oi  02  _  . 

are  defined  as  r0  =  j-j—  ,  r1  =  - —  ,  r  =  - —  ,  n  ■  circular 

41 03  i03  X03 
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Fig.  4 


(Meirovltch,  1970:  286 


.  Strutt  Diagram 
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orbital  frequency,  and  u  *  mass  (trolley)  motion  frequency. 
The  stability  of  the  solution  to  Eq.  (45)  will  therefore 
depend  on  the  station’s  inertia  shape  (r  ,  r2 ) ,  orbit  rate 
to  mass  motion  frequency  ratio  (n/u),  and  the  ratio  of 
effective  trolley  inertia  to  pitch  inertia  rQ.  Since  this 
study  is  only  concerned  with  inertia  parameters ( r  ,  r  )  that 
lie  in  the  classical  Lagrange  region  and  because  a and  qi 
are  a  function  of  r  and  r2,  Chan  (19)  mapped  the  Lagrange 
inertia  region  (Fig.  2)  into  the  Strutt  diagram  (Fig.  4). 

The  result  is  a  polygon  shape  where,  for  a  single  rQ  value, 
the  Lagrange  region  that  is  triangular  in  r  -r  space  maps 
into  a  line  segment  in  the  aT-q^  space  (see  Fig.  5). 

Because  a  point  in  the  two-dimensional  &T-qi  plane  will 
transform  into  a  line  segment  in  the  r^r  plane,  Figs.  6 
and  7  show  how  it  is  easy  to  determine  unstable  space 
station  configurations.  In  Fig.  6,  at  rQ=.25,  points  2  and 
8  represent  what  Chan  (22)  calls  "fence  posts"  enclosing 
points  of  instability.  The  lines  in  Fig.  7  that  correspond 
to  those  points  are  the  boundaries  that  define  the 
approximate  zone  of  instability.  Any  space  station  with  a 
r  -r  combination  that  falls  within  this  zone  will  be 
theoretically  unstable  in  pitch.  In  addition,  as  rQ 
decreases,  the  number  of  points  enclosed  by  the  "fence 
posts"  decreases.  Therefore,  the  width  of  the  instability 
zone  in  the  r  -r  plane  would  also  decrease. 

Each  of  the  diagrams  shown,  have  rQ  values  ranging  from 


STRUTT  DIAGRAM  -  MAPPED  LAGRANGE  REGION  FOR  CU=N) 
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Fig.  5.  Strutt  Diagram  With  Mapped  Lagrange  Region 

For  (w»n) 
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R1  =  (11/13) 


0.0  0.5  1.0  1.5  2.0  2.5  3.0  3. 


_ AT _ 

Fig.  6.  Enclosed  Points  of  Instability 


LRGRRNGE  REGION  WITH  PREDICTED  INSTRBILITY  ZONE 
FOR  (R0=0.25)  RND  (W=N) 


R2  =  (12/13) 


Fig.  7.  Points  of  Instability  Mapped  Into  ri-r*  Space 
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0  to  .25.  The  lower  limit  rQ  =  0  denotes  a  null  trolley 
mass.  The  chosen  upper  limit  rQ=  .25  corresponds  to  the 
elevator  mass  being  equal  to  a  single  sphere  mass  of  a 
dumbbell  configuration  space  station.  In  addition,  each 
diagram  corresponds  to  a  fixed  trolley  frequency  ( u )  where 
n^1  .  The  lower  limit,  u  =  n  (trolley  frequency  =  orbital 
rate),  sets  a  reasonable  lower  bound  on  trolley  speed. 

By  mapping  the  Lagrange  region  onto  a  Strutt  diagram, 
(Chan,  1986)  showed  that  when  o  was  at  the  lower  limit  only 
a  small  part  of  the  first  unstable  region  was  of  concern 
(see  Fig.  8).  But  when  w  is  increased,  even  that  region 
becomes  of  no  concern.  Fig.  9  shows  that  as  w  increases, 
the  mapped  Lagrange  region  moves  towards  the  origin  of  the 
Strutt  diagram  and  away  from  the  first  unstable  region.  As 
this  occurs,  the  "fence  posts"  no  longer  have  any  points  to 
enclose.  It  is  from  this  analysis  that  Chan  (22)  states 
that  if  u  is  greater  than  1.8  times  the  orbital  rate,  then 
all  space  station  configurations  are  stable  in  pitch  when 

A 

mass  motion  is  along  the  yaw  axis  b  . 

The  validity  of  this  argument  comes  into  question, 
however,  when  recalling  the  derivation  in  Appendix  B  where 
the  pitch  equat ion  is  written  as 

z  +  (aT  +  16q1cos2T  +  16q2cos4T 

+  16q3cos6T  +  16q^cos8T]z  =  0  (46) 


When  (Chan,  1986)  used  Mathieu’s  equation  to  determine  pitch 


Fig .  8.  Overlap  of  Mapped  Lagrange  Region  With 
Unstable  Regions  From  Strutt  Diagram 
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stability  of  the  system,  he  essentially  ignored  the  q2  ,  q3 , 
and  q  terms  that  are  present  in  the  Hill’s  equation.  Chan 
(19)  used  the  fact  that  both  the  Hill  and  Mathieu  equations 
are  linear  to  state  that  each  q.  term  can  be  examined 
independently.  As  a  result,  Chan  (21)  makes  a  reasonable 
statement  that  the  q„  and  q^,  terms  can  be  ignored  since  they 
are  not  functions  of  r  or  r2  and  because  they  are  small  (on 

-  4 

the  order  of  *  10  ).  At  this  point,  however,  Chan  (21) 

states  that  since  q2  is  also  a  function  of  the  variables  r  , 
r 1  ,  r  ,  and  (n/w),  a  Strutt  diagram  can  be  developed  in  the 
aT~q2  plane  (see  Fig,  10).  Chan  (21)  goes  on  to  say  that 
due  to  the  location  of  the  mapped  Lagrange  region  to  the 
horizontal  axis,  no  additional  unstable  configurations  are 
added  to  those  obtained  from  the  Strutt  diagram  in  the  aT~q1 
plane.  As  a  result,  Chan  (21)  also  drops  the  q2  term. 

Conducting  the  pitch  analysis  by  using  the  full  pitch 
equation  (all  q.  terms)  demonstrates  that  the  preceding 
arguments  are  not  totally  valid.  The  next  section, 
detailing  Hill’s  method,  shows  how  the  zones  of  instability 


for  specific  r0  and  w  values  in  the  r^r  pi 


ane  are 


different  than  those  predicted  by  Mathieu  Theory. 


Hill’s  Equation 

Determining  the  boundaries  of  the  instability  zones 
using  Hill’s  equation  is  quite  different  from  the  methods 


used  by  (Chan,  1986).  The  main  difference  is  that  a  Strutt 
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Pig.  10.  Strutt  Diagram  In  Av-Q*  Space  Vlth 
Mapped  Lagzange  Region  For  (v=n) 
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diagram  can  no  longer  be  used  to  determine  which  space 
station  configurations  result  in  unbounded  solutions  to  the 
pitch  equation.  The  reason  for  this,  discussed  in  the 
theory  section,  is  that  the  parameter  space  representing 
bounded  and  unbounded  solutions  would  be  in  several 
dimensions.  Fortunately,  Hill’s  equation  is  a  well  studied 
differential  equation  with  solutions  defined  by  (Hayashi, 
1964:341). 

For  the  unbounded  solution  corresponding  to  the  first 
Mathieu  unstable  region,  the  characteristic  exponent  p, 
shown  in  Eq.  (44),  is  given  by 

P  =  6 1  s  i  n  (  2<.T  )  +  e1e2sin^2  o)  +  ^  0203sin(2a) 

"  128  0isin<2a)  +  h  dtd2sin(4o)  +  fgj  e,e3sin(2 o) 


+  ha  ei02sin(2a)  +  hk  e1e3Sin<2o> 

+  h  0203sin(2O)  +  U52  e1e203sin(4a)  + 
in  which  the  parameter  o  is  determined  by 


(47  ) 


K3  +  K1  cos(2c?)  +  K2  cos(4<?) 


(  48a  ) 


where 


K1  *  9,  +  7  V2  *  12  eze,  -  54  +  ife 

-  ih  V2  -  2504  V3  +  h  ele 3 

K2  *  I  e*  ♦  h  e>2  *  Me  ®,V3 

K3  ■  1  -  5  -  l  -  16  03  -  ik  -  ife  VA 

and  where  0Q  =  aT  and  0.  =  8q.  . 

The  task  of  finding  the  instability  zones 


(48b) 


(48c  ) 


(  48d ) 
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if 


initially  appears  to  be  quite  complicated.  However,  the 
boundaries  between  stable  and  unstable  solutions  are  defined 
by  only  two  cases:  (1)  p  =  0,  a  =  0°,  and  (2)  M  =  0,  o  - 

-90°.  Substitution  of  these  o  values  into  Eq.  (48a)  yields 
eQ  =  K3  +  K1  +  K2  (for  o  =  0°)  (49a) 

eo  =  K3  -  K1  +  K2  (for  a  =  -90°)  (49b) 

Since  0Q ,  K1 ,  K2 ,  and  K3  are  a  function  of  aT  and  q.  , 
and  the  values  of  aT  and  qj  are  a  function  of  r  and  r  ,  the 
boundaries  of  the  instability  zone  in  the  r^-r  plane  are 
defined  by  r  ,  r  combinations  that  satisfy  Eqs.  (49a)  and 
(49b)  for  a  given  rQ  and  u. 

Hayashi  (342)  goes  on  to  define  the  unbounded  solutions 
for  other  unstable  regions,  but  they  lie  outside  our  area  of 
interest  and  are  not  included  in  this  analysis.  Therefore, 
a  comparison  now  can  be  made  with  the  results  of  Mathieu’s 
equation  and  to  the  nonlinear  program. 


Results/Observations 


This  section  compares  the  analytical  results  generated 
by  Mathieu’s  and  Hill’s  equations  with  the  numerical  results 
generated  by  the  computer  code  shown  in  Appendix  D.  Each  of 
the  methods  are  compared  at  mass  motion  frequencies  (u) 
ranging  from  1.0*n  to  2.0*n.  In  addition,  at  each  <*>,  the 
methods  will  be  compared  at  two  rQ  values. 

For  comparison  purposes,  the  test  case  run  by  the 
nonlinear  code  is  the  same  case  used  by  (Chan,  1986).  A  low- 
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earth  orbit  of  150  km  was  used  which  yielded  a  circular 

orbit  rate  n  =  .001197  sec  1  and  a  period  of  5250  seconds. 

All  data  points  were  run  for  a  period  of  1000  orbits, 

although  unstable  motion  usually  appeared  by  100  orbits. 

Only  at  points  close  to  the  boundary  between  stable  and 

unstable  motion  when  V  is  small  did  the  simulation  take  more 

than  500  orbits,  as  would  be  expected,  to  go  unstable. 

Examples  of  the  output  generated  by  the  computer  code  for 

some  of  the  data  points  shown  if  Fig.  6  are  shown  in  Figs. 

11  and  12.  Each  figure  shows  the  progression  from  stable 

motion  to  unstable  and  then  back  to  stable  as  r  is  varied 

along  a  constant  r  line.  The  fact  that  is  evident  from 

o 

these  figures  is  that  the  points  predicted  by  the  computer 
code  to  mark  the  beginning  and  end  of  unstable  motion  do  not 
match  the  prediction  by  Mathieu’s  equation.  This 
observation  is  illustrated  more  fully  in  the  comparisons 
discussed  next. 

The  following  eight  cases  detail  the  comparison  of  each 
method  as  u  and  rQ  vary:  (1)  u  =  1.0*n,  (rQ=.25,  rQ=. 175 ) ; 
(2)  o  =  1 . 2*n ,  (r0=.25,  rQ= .175) ;  (3)  «  =  1.4*n,  (r0=.25, 
ro=.20);  (4)  w  =  1.8*n,  (rQ=.25,  rQ=.20).  These  cases  are 
illustrated  in  r  -r  space  in  Figs.  (13)  through  (20). 

Fig.  13  shows  three  predicted  instability  zones  mapped 
into  the  classical  Lagrange  stability  region.  What  is 
noticeable  is  that  the  Hill’s  analytical  prediction  falls 
within  the  Mathieu  predicted  instability  zone.  In 
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orbits 

(1) 


— I  —90 
100  0 


orbits 

(3) 


Pig.  11.  Computer  Runs  (w»n,  ro»0.25) 


Pig.  12.  Computer  Runs  (v=n,  ro*0.175) 
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(11/13) 
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(11/13) 


LAGRANGE  RE6I0N  WITH  PREDICTED  INSTABILITY  ZONE 
FOR  (R0=0.25)  AND  (W=1.2»N) 


Pig.  15.  Instability  Zones  For  (ro*.25)  and  (v*1.2*n) 


LAGRANGE  REGION  WITH  PREDICTED  INSTABILITY  ZONE 
FOR  (R0=0. 175)  AND  (W=1.2*N) 


1.2 

0.4 

M  t.O  1.2 

1.4 

1 .4 

1.t 

2.1 

R2  =  (12/13) 

Fig.  16.  Instability  Zones  For  (ro*.175)  and  (v«1.2*n) 
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(11/13) 


LAGRANGE  REGION  WITH  PREDICTED  INSTABILITY  ZONE 
FOR  (R0=0.25)  AND  (W=1.4*N> 


Fig.  17.  Instability  Zones  For  (ro*.25)  and  (w=1.4*n) 
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(TT7I3)  I  R1  =  (11/13) 


LAGRANGE  REGION  WITH  PREDICTED  INSTABILITY  ZONE 
FOR  (R0=0. 25)  f*C  CW=1 .8-N) 


Fig.  19.  Instability  Zones  For  (ro-.2S)  and  (w=»1.8*n) 


LAGRANGE  REGION  WITH  PREDICTED  INSTABILITY  ZONE 
FOR  (R0=O.2Q)  AND  (W=1.8*N> 


Fig.  20.  Instability  Zones  For  (ros.20)  and  (w»1.8*n) 
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fact,  both  methods  share  a  common  lower  boundary  due  to 


the  fact  that  Mathieu’s  equation  is  a  specialized  form  of 
Hill’s  equation.  Fig.  13  also  reveals  that  both 
analytical  methods  miss  the  same  portion  of  the  instability 
zone  predicted  by  the  computer  code.  It  should  be  noted, 
however,  that  although  Fig.  13  gives  the  impression  that 
the  predictions  by  Hill’s  method  and  the  nonlinear  program 
share  a  common  upper  boundary,  this  is  not  the  case.  Due  to 
the  scale  of  the  graph,  the  difference  between  boundaries  is 
not  obvious.  As  o  increases  or  rQ  decreases,  however,  the 
difference  becomes  more  pronounced. 

Fig.  14  shows  that  the  common  lower  boundary  of  Hill’s 
and  Mathieu’s  predicted  instability  zones  does  not  move  much 
from  its  position  in  Fig.  13  as  rQ  decreases.  All  of  the 
other  boundaries,  however,  will  tend  to  collapse  around  this 
line.  This  tendency  continues  as  rQ  decreases  further 
until  all  that  is  left  is  a  line  at  rQ*  0.  At  rQ=  0  (the 
null  trolley  mass)  no  instability  remains.  It  should  be 
noted  that  this  collapsing  phenomenon  occurs  in  all  the 
cases . 

Another  phenomenon  observed  is  that  as  u  increases,  the 
instability  zones  appear  lower  in  the  Lagrange  region.  This 
downward  migration  of  the  instability  zones  can  be  seen  in 
the  progression  of  graphs  shown  in  Figs.  (13)  through  (20). 
Figs.  (19)  and  (20),  for  example,  show  that  at  w  =  1.8*n  all 
that  is  left  of  the  original  instability  zones  is  the  tip  of 
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the  Lagrange  region.  When  u  >  1.8*n,  r.o  instability  zones 
rema i n . 

The  predicted  instability  zones  do  not  all  disappear 
at  the  same  time.  The  zone  predicted  by  Hill’s  equation, 
for  instance,  disappears  by  the  time  u  =  1.3*n.  As  a 
result,  Figs.  (17)  and  (18)  show  only  two  zones  remaining. 
When  w>  r  1.7*n,  the  zone  predicted  by  the  nonlinear 
equations  vanishes  leaving  Figs.  (19)  and  (20)  with  only  one 
zone  . 

Based  on  this  information,  neither  analytical  method 
provides  an  excellent  prediction  of  what  the  nonlinear 
equations  are  saying.  Depending  on  the  purpose  of  their 
use,  however,  the  analytical  methods  maybe  adequate  for 
determining  zones  of  pitch  instability  in  the  Lagrange 
region.  Therefore,  a  summary  of  the  results/observations  is 
given  below. 

For  1.0*n  S  «  <  1.3*n,  both  analytical  methods  predict 
and  miss  the  same  portions  of  the  nonlinear  instability 
zone.  In  addition,  Hill’s  method  fails  to  predict  any  new 
instability  zones  that  have  not  already  been  predicted  by 
Mathieu’s  method.  The  conservative  nature  of  Mathieu’s 
method,  however,  eliminates  an  additional  group  of  space 
station  configurations  that  Hill’s  method  does  not. 

For  1.3*n  s  <j  <  1.7*n,  Mathieu’s  method  predicts  a  much 
larger  instability  zone  than  the  one  predicted  by  the 
nonlinear  equations.  However,  the  nonlinear  zone  hovered 
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around  the  lower  boundary  of  Mathieu’s  zone.  Given  that 
Hill’s  method  disappeared,  Mathieu’s  analytical  solution  is 
better  than  nothing  because  it  predicts  the  approximate 
neighborhood  of  the  nonlinear  zone. 

For  1.7*n  £  m  <  1.9*n,  Mathieu’s  method  still  predicts 
an  instability  zone  when,  according  to  the  nonlinear  program, 
an  instability  zone  no  longer  exists. 

And  finally,  for  u  £  1.9*n,  no  instability  zones  are 
predicted  by  any  method. 


VI .  Roll /Yaw  Analysis 


The  objective  of  this  chapter  is  to  obtain  an  analytical 
solution  to  the  coupled  roll/yaw  equations  using  the  Method 
of  Multiple  Scales.  This  method  is  used  because  the 
analytical  technique  in  Chapter  V  is  not  applicable  to 
two-dimensional  problems.  The  results  are  compared  to  the 
numerical  solutions  of  the  nonlinear  system  for 
verification. 

The  chapter  begins  by  reviewing  the  theory  behind  using 
perturbation  techniques  to  solve  parametrically  excited 
systems.  The  Method  of  Multiple  Scales  is  then  introduced, 
followed  by  a  discussion  on  the  existing  resonance 
conditions.  The  remainder  of  the  chapter  discusses  the 
results  obtained  and  the  comparison  to  the  numerical 
solutions . 


Theory 


The  linearized  roll  and  yaw  equations  of  attitude  motion 
were  shown  to  be  coupled  in  Chapter  IV.  Re-writing  Eqs .  (32) 
and  (33)  yields 
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(  50  ) 


0(51) 


where  each  coefficient  containing  I  is  a  function  of  time 


due  to  mass  motion.  Since  the  mass  motion  is  periodic  and 
occurs  in  an  environment  where  the  gravity-gradient  torque 
is  the  dominant  factor  affecting  the  attitude  stability  of 
the  space  station,  the  equations  governing  attitude  motion 
are  parametrically  excited.  The  presence  of  this  small 
disturbance  term  has  a  non-negl igible  cumulative  effect 
because  it  acts  over  a  long  time. 

There  are  several  techniques  available  to  provide  a 
systematic  way  of  estimating  this  cumulative  effect.  These 
techniques  are  generally  referred  to  as  "perturbation 
methods"  because  the  nonlinear  equations  representing  the 
non  autonomous  dynamical  system  are  "quasi-linear” .  This 
means  that  the  equations  can  be  separated  into  one  part 
containing  linear  terms  and  a  second  part  containing  small 
nonlinear  terms  known  as  perturbations  (Meirovitch,  1970: 
294).  The  perturbation  technique  used  in  this  study  is 
known  as  the  Method  of  Multiple  Scales  or  Two-Timing. 

Method  o f  Multiple  Scales 

Appendix  C  details  the  expansion  which  shows  Eqs.(50) 
and  (51)  are  of  the  form 

Uj+  1 a1ut  =  0  (52) 

u2  -^Ujt  a2u2+  2e(f2iui+  f22u2  >sin(  2wt  )  (53) 

+  2e(f23U1+  f2,u2)cos(2ut)  =  0 
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where 
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w* 

3 

1 _ 

y  1 

L“«J 
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(  54  ) 


and 


w  ■  frequency  of  mass  motion 
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When  £  =  0,  we  have  a  system  of  linear  homogeneous 
differential  equations  with  known  solutions.  The  role  of 
the  parameter  £,  therefore,  is  to  separate  out  the  small 
terms  in  Eq.(53).  Meirovitch  (295)  states  that  the  presence 
of  such  a  parameter  associated  with  the  small  terms  enables 
us  to  effect  the  transition  between  the  known  and  desired 
solutions . 

At  this  point,  the  fundamental  perturbation  technique 
assumes  a  solution  to  Eqs.  (52)  and  (53)  in  the  form  of  a 


46 


power  series  in  c  given  by 


u  =  u  +  £u  +  £  u  +  •••  .  (56) 

m  m  0  m  1  m  2 

However,  the  method  of  multiple  scales  introduces  new 
independent  variables  according  to 

V  ^  for  n  =  o,  1,  2,  •••  (57) 

where  T0=  t  ;  and  Tt=  et  . 

The  underlying  idea,  according  to  (Nayfeh,  1979:  56),  is  to 
arrive  at  a  more  convenient  method  of  dealing  with  certain 
dynamical  systems  by  considering  "the  expansion  representing 
the  response  to  be  a  function  of  multiple  independent 
variables,  or  scales,  instead  of  a  single  variable."  As  a 
result,  the  assumed  solution  becomes 

u  =  u  ( T  ,  T  )  +  Cu  ( T  ,  T  )  +  •••  (58) 

m  mOOl  m  1  0  1 

where  the  independent  time  scales  are  limited  to  TQand  T 
because  a  first  order  expansion  is  adequate  for  this  study. 


Expanding  the  time  derivatives  of  Eq .  (58)  yields 


0  nO 

where  the  operators 


u  =  Du  +  eD  u  +  eD  u  +  e  Du 
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Similarly, 


u  =  D^u  +  2£D„D,u  +  eDfu  ,+  e?Dfu  e^D^D.u  + 

m  0  » 0  0  1  m  0  Oml  ImO  01ml 


•  (  61  ) 

Substituting  Eqs.  (53)  and  (61)  into  Eqs.  (52)  and  (53)  and 

equating  coefficients  of  like  powers  of  £  yields 
o 

Order  £ 


2 

Du  +  XDu  +  au  = 

010  1020  110 

0. 

(  62a  ) 

Du  -  X  D  u  +  a  u  = 

0  20  20  10  2  20 

0  . 

(62b) 

Order  e 1 

2 
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(  63a  ) 

2 

Du-XDu  +  a  u  =  -2D  Du  -  X  D  u 
0  21  2011  2  21  01  20  2110 


■2<f21U10+  f22D0U20)Sin(2tJt) 


(  63b  ) 


—  2 ( f  Du  +  f  u  ) cos ( 2ut ) 

23  0  10  24  20;  '  ' 

According  to  (Nayfeh,  1979:  332)  the  solutions  to  Eqs.  (62a) 
and  (62b)  can  be  written  in  the  form 
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where  and  u are  the  solutions  to  the  order  £ 


characteristic  equation 


4  2 

w-(a+a+XX)o+aa  =0. 

l  2  1  2  ’  12 


(  64b) 


(65  ) 


It  is  assumed  that  ^  and  w2  are  real  and  positive,  and  that 


w 


Now  that  the  solutions  to  the  linear  homogeneous 
differential  Eqs.  (62a)  and  (62b)  are  known,  we  can 
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substitute  them  into  Eqs.  (63a)  and  (63b)  to  get 
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where  A1  and  A2  are  derivatives  of  A^  and  A2  with  respect  to 
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where  A  and  A2  are  the  complex  conjugates  of  A ^  and  A2 . 


Resonance 

Fortunately,  it  is  not  necessary  to  find  the  solutions 
to  Eqs .  (66)  and  (68)  in  order  to  determine  the  system’s 

stability.  This  study  needs  only  to  look  at  the  terms 
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.. 

from  Eq .  (68).  The  reason  for  this  stems  from  Eq.  (65)  and 

the  determination  of  (*>1  and  <J2 .  The  frequencies  and  i*>2 
are  the  resonant  frequencies  of  the  unperturbed  equations  of 
motion.  When  the  mass  motion  frequency  w  takes  on  a  value 
such  that  one  of  the  exponents  in  Eq.  (66)  goes  to  a 
resonant  frequency  (i.e.  iu1T(J  or  iu2T0 ) »  secu^ar  terms  will 
appear  in  the  solution.  Since  secular  terms  grow 
indefinitely  with  time,  large  oscillations  in  attitude  or 
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totally  unstable  motion  can  occur.  The  resonance  phenomenon 
in  this  study  occurs  when  a  mass  on  the  space  station  moves 

A 

along  the  b1  axis  with  a  frequency  (u)  that  meets  one  of  the 
following  conditions:  (a)  u  *  ,5(U2+  )  !  (b)  u  w  .5(W2~ 

)  ;  (  c  )  u  ;  or  (  d  )  u  w  .  For  example  ,  when  w  is  at 


T  i  (  )  Tq'|  r  i 

condition  (a),  [e  J  becomes  [e 

‘U2Tol 

which  equals  ^e  J  which  is  a  resonai 


(  2* ( . 5u  +u ) -w  ) 


which  equals  ^e  J  which  is  a  resonant  frequency. 

This  problem  of  resonance  and  gravity-gradient  stability 
has  been  studied  before  by  (Breakwell  and  Pringle,  1965). 

In  their  study,  Breakwell  and  Pringle  (307)  identified  orbit 
eccentricity  instead  of  mass  motion  as  the  forcing  function. 
Fig.  (21)  contains  Hotted  lines  which  represent  what  they 
call  "external  resonance"  loci  in  the  i'1,r2  inertia 
parameter  space.  The  authors  stated  that  if  a  spacecraft 
had  a  configuration  that  fell  on  an  external  resonance  line, 
orbit  eccentricity  would  eventually  cause  unbounded  attitude 
motion . 

The  types  of  resonance  lines  described  by  (Breakwell  and 
Pringle,  1965)  also  appear  in  this  study.  The  reason  for 
this  is  that  <«>1  and  <*>2  are  a  function  of  r  Q,  r^  and  r2 .  As 
a  result,  each  resonance  condition  translates  into  constant 
value  lines  (loci)  in  the  Lagrange  region.  To  illustrate 
this  concept,  Figs.  (22)  through  (25)  show  the  constant 
value  lines  created  by  the  resonance  conditions  (a)  through 
(d)  at  rQ=  .25.  Figs.  (26)  and  (27)  show  how  constant  value 
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(Break we 11  6  Pringle,  1965:  307) 


Pig.  21.  Internal  and  External  Resonance  Lines  in 
ri,ra  Inertia  Parameter  Space 


LAGRAN6E  STABILITY  REGION  -  RESONANCE  LINES  AT  <RO=.25> 
FOR  .5»(W2+W1)  =  MASS  MOTION  FREQ.  (W) 


Fig.  22.  Resonance  Lines  For  .5(W2+W1) 
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LAGRflNGE  STABILITY  REGION  -  RESONRNCE  LINES  ftT  (RO=.25) 
FOR  M2  =  MRSS  MOTION  FREQ.  (W) 
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Fig.  24.  Resonance  Lines  For  W2 
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Fig.  26.  Resonance  Lines  For  .5(W2+W1)  As  A  Function  of  ro 
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lines  for  conditions  (a)  and  (c)  vary  with  rQ. 

The  lines  in  each  figure  represent  where  in  the  Lagrange 
region  that  the  designated  function  equals  a  constant  value. 
The  values  are  varied  in  each  figure  to  show  the  progression 
of  the  lines.  For  example,  when  a  mass  moves  on  a  space 
station  with  a  frequency  (u)  equal  to  1.0*n,  the  constant 
value  line  for  ,5*(u  +w  )  in  Fig.  (22)  that  corresponds  to 
1.0*n  becomes  a  resonance  line.  In  addition,  as  u  is  varied 
so  is  the  location  of  the  corresponding  line.  It  should  be 
noted,  however,  that  only  two  of  the  four  resonance 
conditions  listed  previously  result  in  resonance  lines  of 
concern.  Figs.  (23)  and  (25),  for  example,  show  the 
constant  value  lines  for  the  resonance  conditions  .  5*  ( ) 
and  u  .  The  figures  reveal  that  in  each  case  the  functions 
do  not  exist  in  the  Lagrange  region  for  values  of  w  greater 
than  1.0*n.  Figs.  (22)  and  (24),  however,  reveal  constant 
value  lines  for  values  of  u  between  1.0*n  and  2.3*n  for  the 
conditions  (a)  u  *.  5*  ( U)2+(J1  >  and  (c)  w  «  u>2 ,  Since  the 
lower  limit  of  u  was  determined  in  Chapter  V  to  be  1.0*n, 
this  study  will  only  look  at  the  resonance  lines  generated 
by  conditions  (a)  and  (c). 

In  addition  to  the  downward  progression  of  the 
resonance  lines  due  to  increasing  w,  Figs. (26)  and  (27)  show 
a  downward  shift  due  to  a  decrease  in  rQ.  The  figures  show 
that  the  shift  is  more  noticeable  in  condition  (c)  and  when 
w  is  high. 
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In  order  to  verify  that  resonance  conditions  (a)  and  (c) 


indeed  lead  to  unstable  motions,  computer  runs  were 
conducted  at  a  number  of  specific  points  on  the  different 
constant  value  lines  for  each  condition.  Two  points,  one 
from  each  condition,  are  shown  in  Figs.  (28)  and  (29).  The 
points  illustrated  were  picked  so  that  they  would  lie 
outside  the  pitch  instability  zones  described  in  the 
previous  chapter,  in  order  to  assure  that  the  results  were 
directly  attributed  to  the  resonance  phenomenon.  Fig.  (28) 
shows  that  the  pitch,  roll,  and  yaw  attitude  motion  for  a 
point  on  the  1.0*n  line  of  .Sfcfc^+o^)  is  unstable.  Fig. 

(29)  reveals  the  same  results  for  a  point  on  the  1.4*n  line 
of  u2. 

It  should  be  realized,  however,  that  the  unstable  motion 
associated  with  these  resonance  lines  is  not  confined 
specifically  to  the  lines  themselves.  As  discussed  earlier, 
the  perturbation  parameter  £  associated  with  the  small 
nonlinear  terms  enables  us  to  effect  the  transition  between 
the  known  and  desired  solutions.  This  transition  acts  as  a 
band  around  the  resonance  lines  to  mark  the  change  from 
stable  to  unstable  motion.  To  illustrate  this  fact,  Fig. 

(30)  enlarges  the  area  around  the  point  used  for  Fig.  (28). 
Two  points  are  chosen  at  varied  distances  from  the  line.  A 
computer  run  is  conducted  using  the  nonlinear  program  to  see 
what  happens  to  the  stability  of  the  equations  at  each 
point.  Figs.  (31)  and  (32)  reveal  that  point  (1)  is  still 
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Fig.  28. 

Computer  Run  For  Resonance 
Line  . 5(W2+W1)  at  (v=1.0*n) 
(ri.=.6841#ra=.76) 


Fig.  29. 

Computer  Run  For  Resonance 
Line  W2  at  (w=1.4*n) 
(ri=.73,r*=.85) 


58 


(11/13) 


Fig.  30. 


Sample  Bandwidth  For  Resonance  Line 
. 5 ( V2+V1 )  at  ( v=l . 0*n ) 
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Fig.  31. 

Fig.  32. 

Computer  Run  For 

Computer  Run  For 

Bandwidth  Pt(l) 

Bandwidth  Pt(2) 

raw 


within  the  transition  band  where  as  point  (2)  is  not. 
Although  both  points  were  run  for  1000  orbits,  only  a 
portion  of  the  data  is  presented  so  that  the  trends  can  be 
observed . 

Mathematical  approximation  of  the  bandwidth  around  the 
resonance  lines  is  possible  and  procedures  are  laid  out  in 
(Nayfeh,  1979).  However,  the  bandwidths  are  a  function  of  e 
and  therefore  depend  on  rQ  and  v 2 .  Although  rQ  is  a  chosen 
constant,  the  value  of  r2  varies  along  the  resonance  lines. 
Therefore,  the  bandwidth  will  have  a  tapered  appearance  with 
the  width  depending  upon  the  location  in  the  Lagrange 
region.  For  this  study,  it  is  considered  sufficient  to  say 
that  a  band  exists  around  the  resonance  lines  and  that  it 
should  be  considered  when  making  conclusions.  Determination 
of  the  actual  bandwidth  is  left  for  future  studies. 

Results 

The  main  idea  that  this  chapter  attempted  to  present  was 
that  in  the  Lagrange  region,  resonance  lines  exist  which 
represent  space  station  configurations  that  result  in 
unstable  attitude  motion  for  specific  mass  motion 
frequencies  (u).  The  lines  that  would  cause  concern  at 
specific  u’s  are  presented  in  Figs.  (33)  through  (36).  It 
should  be  noted  that  as  w  increases,  the  resonance  lines 
migrate  towards  the  bottom  of  the  Lagrange  region.  By  the 
time  u=2.3*n,  no  resonance  lines  remain.  If  a  space  station 


Fig.  33.  Resonance  Lines  When  (w=1.0*n) 
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LAGRANGE  STABILITY  REGION  -  RESONANCE  Llt€S  AT  <A0=.25> 


HASS  NOTION  FREQ.  (W>  =  1.8*»<n> 


Fig.  35.  Resonance  Lines  When  (w=1.8*n) 


63 


A 

is  designed  such  that  mass  motion  will  occur  along  the  b1 
axis  with  a  1.0*n  £  u  £  2.3*n,  then  configurations  which  lie 
in  the  neighborhood  of  the  resonance  lines  should  be  avoided 
or  at  least  analyzed  further  with  appropriate  consideration 
given  towards  compensation. 
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In  order  to  make  conclusions  and  recommendations  on  the 
analyses  provided  in  Chapters  V  and  VI,  a  composite  of  the 
results  obtained  in  both  chapters  is  formed.  To  accomplish 
this  task,  some  ground  rules  are  established. 

In  the  pitch  analysis  conducted  in  Chapter  V,  two 
analytical  methods  and  the  numerical  solutions  were 
compared.  Each  method  predicted  a  different  instability 
zone  and  each  method  existed  throughout  a  different  range  of 
mass  motion  frequencies  (w).  Since  neither  analytical 
method  exhibited  superior  qualities  relative  to  the  other, 
no  preferred  method  was  picked.  Therefore,  in  order  to 
simplify  the  overall  analysis,  the  Mathieu  analytical  method 
is  used  to  represent  the  pitch  analysis  because  it  was  the 
most  conservative  and  it  existed  over  the  largest  range  of 
applicable  u  values. 

In  the  roll/yaw  analysis  conducted  in  Chapter  VI, 
resonance  lines  were  determined  for  different  values  of  w. 

In  addition,  a  band  around  each  resonance  line  was  found  to 
exist.  Because  no  analytical  method  is  developed  in  Chapter 
V  to  determine  a  specific  bandwidth,  the  composite  analysis 
uses  only  the  actual  resonance  lines. 

When  the  findings  from  Chapters  V  and  VI  are  combined 
under  the  ground  rules  stated  above,  the  result  is  a 


composite  of  instability  regions  as  a  function  of  u.  Figs 
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# 


(37)  and  (38)  show  all  the  r),r2  combinations  that  lead  to 
unstable  motion  for  the  given  mass  motion  frequency  (u). 

Conclusions 

A  detailed  discussion  of  the  different  analytical 
methods  used  in  the  pitch  and  roll/yaw  analyses,  along  with 
the  presentation  of  results,  is  located  at  the  end  of  each 
of  the  respective  chapters.  Overall,  the  following 
conclusions  can  be  drawn: 

1.  For  the  problem  of  a  simple,  gravity-gradient  stabilized 
rigid  body,  the  addition  of  mass  motion  results  in  time 
varying  moments  of  inertia  being  introduced  into  the 
equations  of  motion. 

2.  The  mass  motion  must  be  periodic  and  continuous  over 
many  orbital  periods  to  induce  parametric  excitation. 

3.  A  composite  of  the  pitch  and  roll/yaw  instability 

A 

predictions  when  the  mass  motion  is  along  the  b^  axis  yields 
the  following:  (1)  for  1.0*n  s  w  <  1.5*n,  there  exists  a 

zone  of  pitch  instability  along  with  two  roll/yaw  resonance 
lines;  (2)  for  1.5*n  i  u  i  1.8*n,  there  exists  a  zone  of 
pitch  instability  along  with  one  roll/yaw  resonance  line; 

(3)  for  1.8*n  <  w  <  2.3*n,  there  exists  one  roll/yaw 
resonance  line;  (4)  for  £  2,3*n,  no  predicted 
instabilities  exist. 

4.  Mathieu  and  Hill  analyses  of  the  linearized  pitch 
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Fig.  37.  Composite  of  All  Unstable  Space  Station 
Configurations  at  (w=1.0*n)  and  (ro=.25) 


Fig.  38.  Composite  of  All  Unstable  Space  Station 
Configurations  at  (w=1.4*n)  and  (ro=.25) 
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equation  provide  a  1st  order  prediction  of  unstable  space 
station  configurations  resulting  from  mass  motion.  These 
methods  should  be  used  as  a  guide  to  establish  a  more 
precise  location  of  unstable  configurations  using  the 
nonlinear  equations. 

Re commend at  ions 

The  following  recommendations  for  further  study  are 
suggested : 

1.  Study  of  the  problem  of  mass  translation  along  other 

A  A 

axes  (i.e.  b^  or  b3  )  may  pro%'ide  the  potential  for 
additional  information  gain. 

2.  Introduction  of  other  forcing  functions,  like 
eccentricity  or  the  rotation  of  internal  mass,  to  the  rigid 
body  problem  is  recommended. 

3.  Study  of  a  more  realistic  mass  motion  problem  that 
includes  flexibility  effects  along  with  eccentricity  for  a 
specific  configuration  (i.e.  proposed)  is  recommended. 

4.  Further  study  is  recommended  to  better  predict  the 
roll/yaw  instability  regions  by  determining  a  mathematical 
approximation  to  the  bandwidths  around  the  resonance  lines. 
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Appendix  A 


Equat ions  o f  Mot  ion  Using  Eul e r  Parameters 


As  stated  in  Chapter  III,  quaternions  are  a  more 
efficient  parameterization  of  the  attitude  position  of  a 
spacecraft  than  Euler  angles.  The  reason  for  this  is  that 
they  are  more  compact  than  the  direction  cosine  matrix. 
Quaternions  require  less  computational  time  in  simulations 
since  there  are  only  four  parameters  instead  of  nine  and 
because  they  involve  fewer  trigonometric  functions. 

(Wertz,  1980:  512)  describes  the  time  differentiation  of 
quaternions  as 

=  <l/2)Za  ( A- 1  ) 

where  Z  is  the  skew-symmetric  matrix 
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W 

expressed 


the  components  of  the  angular 
in  the  principal  body  basis.  Since, 


IB  10  OB 
W  =  W  +  K 


(  A-4  ) 


A  n  +  w 

1  J  u 

A  n  +  w 

2  J  v 

A  n  +  w 

3  3  w 


(  A-  5  ) 


Solve  for  w  ,  w  ,  and  w  ,  and  substitute  them  into  Eq 
u  v  w 


( A- 1 )  to  get 


qt  =  (l/2)(q2Ww-q3Wv+q4Wu) 

q2  =  (l/2)(-qiWw+q3Wu+q<Wv) 

q3  =  < l/2)(qiWv-q2Wu+q4Ww) 

S  =  (l/2)*-qiWu-q2Wv-q3Ww) 


( A-6a  ) 


( A-6b ) 


( A-6c  ) 


( A-6d  ) 


Eqs .  ( A—  6 )  can  be  rewritten  as 

q,  =  ( l/2)[q2(w3-A33n)-q2(w2-A23n)+q4(wi-A13n) ]  (A-7a) 

q2  =  { 1/2  )  f -qj ( w3-A33n)+q3( wi-A1 3n )+q4 ( w2-A23n ) ]  (A- 7b) 


q3  =  (l/2)[qi(w2-A23n)-q2(wi-A13n)+q4(w3-A33n)] 


( A- 7c  ) 


q  ^  =  (l/2)[-q1(w1-A13n)-q2(w2-A23n)-q3(w3-A33n)]  (A-7d) 

Eqs.  (A-7)  are  the  equivalent  set  of  kinematic  equations  to 
be  solved  simultaneously  with  Eqs.  ( 1 ) »  (2),  and  (3)  instead 

of  Eqs.  (17).  In  addition  to  Eqs.  (A-7),  the  following 


relations  are  needed 


q,  =  (1/4q<)(A23-A32) 

q  2  =  (1/4q<><A3l-A13) 

q  3  =  d/4q4)(A12-A21) 


( A-8a) 
( A-8b) 


( A-8c ) 
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(  A- 8d  ) 


q4  =  (1/2)(1+A1i+A22+A33+A,4) 


.  5 


which  give  Euler  parameters  in  terms  of  the  elements  of  the 
direction  cosine  matrix  A.  Similarly,  A  in  terms  of  q  can 
be  seen  to  be 


A  (  y  !  = 


q1-q2-q3+q4 

2 ( q iq2-q3q4  ) 

2(qiq3  +  q2q<  > 


2(q,q2+q3q4 ) 
2  2  2  2 

"qi+q2"q3+q4 

2(q2q3-qiq4 ) 


2(qiq3"q2q4) 

2(q2q3+qiq4> 

2  2  2  2 

-qi'q2+q3+q4  - 


(  A-9  ) 
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ix  B 


Expansion  of  Periodic  Coefficient 


The  moments  of  inertia  about  the  composite  center  of 
mass  of  a  body  that  is  symmetric  about  the  yaw  axis  and  has 
mass  movement  along  the  yaw  axis,  can  be  described  by 


*1  =  *01 

12  =  I02  +  M  [  (L/2  )sin(ot  )  ]: 

13  =  I03  +  M  [ ( L/2 >sin(wt ) V 


( B- 1  a ) 


( B-lb  ) 


( B-lc  ) 


where  M  is  the  reduced  mass  of  the  moving  element.  If  the 
space  station  assumes  the  form  of  a  dumbbell  satellite,  for 
example,  M  =  [mm  /(m  +m  )]  where  m  is  the  mass  of  the 

e  8  e  s  e 

moving  element  and  ms  is  the  mass  of  one  of  the  spheres. 

!01,  I02>  anc*  I  o  3  1  are  the  principal  moments  of  inertia 

about  the  center  of  mass  of  a  generic  space  station  without 
the  moving  mass.  Time  derivatives  of  Eqs.  ( B— 1 )  yields 

I  =  i2=  i3=  [  wML2sin  (  ut  )  cos  (  ut )  ]  (B-2) 

Inserting  Eq.  (B-2)  into  Eq.  (37)  gives  us 


P( t  )  =  .5 


.5u)ML2sin(o)t)cos(cJt) 
IQ3  +  .  25ML2sin2 ( wt ) 


(B-3  ) 


Squaring  P(t)  gives 


p2(t)  _  25  •  2  5a2M2L*sin2  ( cot  )cos2  ( <>>t  ) 

L<I03  +  M  [  (  L/2 ) sin( wt ) ] 2  ) 2 


(  B-4 


Taking  a  time  derivative  of  P(t)  yields 
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Pit)  =  .25  ( uML2  )  [ wCcos 2 ( ut ) -wCs in2 ( wt ) 


- . 5^ML2C2sin2 ( ut )cos2 ( wt ) ]  (B-5) 


where  C  =  [I  +  .  25ML2  s  in2  ( (ot  )  ]  *. 

Inserting  Eqs .  ( B— 1 )  into  Eq.  (38)  will  give 


R2(t)  =  3n2C  [ IQ2+  . 25ML 2s in2 ( wt )  -  I  ] 


(  B-6  ) 


Substituting  Eqs.  ( B-4  )  ,  (B-5),  and  (B-6)  into  the 

transformed  pitch  equation,Eq.  (41),  yields 


z  +  (3n  C 


I02+  . 2 5ML2 s in2 ( ut )  -  I 


C2u2M^L4  s  in2  (  ut  )  cos  2  (  ut  ) 
1  b 


]- 
-  \  WML2  [ 


Cwcos  ( wt )  -  (  B- 7 ) 


Cwsin2(ut)  -  ^  wML2C2sin2 ( Wt )cos2 ( wt )  >  z  =  0. 


2  2  - 1 

The  term  C  =  [IQ3+  -25ML  sin  ( «t ) ]  can  be  written  in  the 


form  C  = 


'  1  ' 


1  +  x  ) 


- 1 


(  B-8  ) 


o  y 


2  ML 

where  x  =  rQsin  (wt)  and  rQ  =  -j-r—  • 

'  1 03 


-  1 


By  using  the  binomial  series  for  (1  +  x)  and  only  keeping 
the  first  three  terms  we  get 


=  (d 


(1  -  X  +  X  ) 


2(  1  -  2x  +  3x2  ) 


( B-9  ) 


(B-10) 


The  three  term  approximation  for  the  binomial  series  comes 


from  the  fact  that  rQ  has  a  maximum  value  of  .25,  as 


discussed  on  page  27  of  Chapter  V.  Substituting  (B-9)  and 


( B  —  1 0  )  into  (B-7)  yields 


tel 


IQ2+  . 25ML2sin? ( ut ) 


!01  (1  -  x  +  x  ) 


2  2  4 

—  sin2  {  wt  )cos2  (ut  )  (  1 

16I2 

03 


-  2x  +  3x  )  - 


(B-ll  ) 


2  2  r 

cos2(ut)(  1  -  x  +  x2)  -  sin2{wt)(l  -  x  +  x2  )  - 

1  o  3 


|y —  sin2(wt)cos‘1(<Jt)(l  -  2x  +  3x2)  l  z  =  0. 
1  03  •*  J 


By  defining  the  terms 

,  2  2  4 

u  M  L  2  2 

A  =  - - —  sin  (ut)cos  (wt) 

16I03 

=  u2r2sin2 ( ut ) cos2 ( wt ) 

B'=  flat!  *  „v 

4 1  „  „  0 

0  3 

'  2  2 

C  =  (cos  ( wt )  -  sin  (ut))  =  cos(2ut) 

,  2 

D  =  Tj-y —  sin2  (ut  )cos2  ( wt )  =  2r0sinZ(wt  )cos2(wt ) 

X03 


B  C  =  <>>  rQcos(2wt) 


2  2  2  2 
B  D  =  2u  rQsin  («t)cos  ( ut ) 

Eq.  (B-ll)  becomes 


tel 


Iq2+  • 25ML2sin2 ( ut )  -  I  (1  -  x  +  x2 )  - 


(B-12) 


A  (l-2x  +  3x2)  -  B  C  ( 1-x+x2  )  +  B  D  ( l-2x  +  3x2  )  >  z  =  0. 
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But  3n 


fl  1 

fl  1 

02 

w 

=  3n2r2 

> 

3n2 

01 

=  3n  r 


2  Lj  2  2  2 

and  3n  77 —  sin  ( wt  )  =  3n  r.sin  ( wt ) 

4  1  „  U 

03 

Therefore,  Eq.  ( B— 12)  becomes 

z  +  ^3n2r?  +  3n2rQsin2(ut)  -  3n2r^  -  3n2r2x  -  3n2rQsin2(wt)x 

+  3n2rJx  +  3n2r2-\2  +  3n2rosin2(wt)x2  -  3n2rix2  -  A  +  2A  x  - 

<  '  '  '  '  '  '  2  '  '  '  '  '  '  2I 

3 A  x  -  B  C  +  B  C  x  -  B  C  x  +  B  D  -  2B  D  x  +  3B  D  x  z  =0 

(B-13 ) 

2 

By  recalling  that  x  =  rQsin  ( wt )  and  grouping  like  terms, 

Eq .  (B-13)  becomes 

I  2  2(2  2  2 

z  +  <  3n  r2  -  3n  rt+  3n  rQ  -  3n  r2rQ  +  3n  r,r„|sin  (wt) 


i  +  [3nZr0  -  3n2r2r0  +  3n2riro]sin? 

2)  .  ,  .  2  3  .  6, 

ir0Jsm  ( wt  )  +  3n  rQsin  (w 

(  2  2  2  2I  2  2 

l-w  rQ  +  2w  rQ  sin  («t)cos  (wt) 

I  2  3  1  ^  2  \2  A  1  6  2 

-2w  r0jsin  (wt)cos  (wt)  +  3w  rQ  sin  (ut)cos  (wt) 
|w2r2j sin2 ( wt )cos ( 2wt )  -  ^w2r2j sin* (wt )cos( 2wt ) 


+  |-3n2ro  +  3n2r2r2  -  3n2r 


(B-14  ) 


2  1 

-  w  r0cos(2wt)>  z  =  0. 

By  using  the  trigonometric  identities 
sin2(wt)  =  1/2  -  ( 1/2  Icosl 2wt ) 


sin  (wt)  =  3/8  -  ( 1 /2 ) cos ( 2wt )  +  ( 1/8 )cos( 4wt ) 
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sin  (  wt  )  =  (5/16)  -  (  15/32  )cos(  2wt  )  +  <  3/1  6  )cos  (  4<Jt  )  - 

( 1/32  )cos( 6wt ) 

2  2 

cos  ( wt )  -sin  ( «t )  =  cos(2w>t) 
cos2(ut )sin2(ut  )  =  1/8  -  ( 1 /8 ) cos ( 4«t ) 

and  grouping  like  terms  of  cos(2ut),  cos( 4wt ) ,  cos(6w>t),  and 
cos(  8<Jt  )  ,  Eq  .  (  B  —  1  4  )  can  be  written  as 


+  16qicos(2T)  +  16q2cos(4T)  +  16q3cos(6T) 


+  1  6q  cos ( 8T  )  z  -  0.  ( B- 1 5  ) 

4 


in  which  the  independent  variable  has  been  changed  from  t  to 
T  through  the  definition  T  =  wt  ,  and  the  variables  a^  ,  q^ 
q^ ,  and  q4  are  given  by 

aT  =  (  9/8  )  (  n/u  )  2r^  (  r2-  r1~  1)  -  (  3/2  )  (  n/t*>)  2rQ  (  r2~  r^  1) 

+  3(n/w)2(r2-  ri  )  +  ( 1 5 / 1 6  )  ( n/w ) 2 r*  +  (15/128)rJ 


+  (  1/8  )r30  -  (  1/8  )r2 


( B-16 ) 


q  1  =  (  1/16  )  j^( -3/2  )(  n/u)  2r2(  r2-ri-l  )  +  (  3  /  2  )  (  n/<J  )  2r  Q  (  r2~  r  t  -  1  ) 
-(45/32)(n/u)2r^  -  (3/32)r*  -  (3/8)r*  +  (l/2)r2  -  rQj  (B-17) 
q2  =  ( 1/16) [(3/8)(n/u)2r2(r2-r1-l )  +  ( 9/16  )  ( n/u) 2r^ 


( 3/32  )  r^  +  ( 3/8  )  -  (3/8)rj 


( B-18 ) 


q3  =  (1/16)  -(3/32)(n/w)2r^  -  (3/32)r*  -  (l/8)r^ 


(B-19 ) 


q,  =  (1/16)  [-(3/128)  r^] 


(  B  —  2  0  ) 
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Expansion  of  the  Roll /Yaw  Periodic  Coefficient 


The  roll  and  yaw  equations  of  attitude  motion  were  shown 


in  Chapter  IV  to  be 


I-I-I  I-I 

3  21-  3  2  2 

y  +  - - -  nr  +  - - -  n  y  =  0 

1 1  1 1 


(C-l  ) 


1 3 —  I2-  I,  .  Mg"  2  1  2 1  • 

r  -  - - -  ny  +  4  - - -  n  r  +  — r—  (  r  +  yn  )  =  0 

i2  i2  12 


(  C-2  ) 


where  each  coefficient  containing  I  is  a  function  of  time 
and  contributes  to  exciting  the  attitude  motion  of  the  space 
station.  These  terms  can  be  expanded  in  a  manner  similar 
to  the  pitch  periodic  coefficient  in  Appendix  B. 


Coef  f ic ient 


By  recalling  from  Appendix  B  that 


1 2  =  ! 02  +  M  ( (L/2 )sin(ut ) ] 


I  =  i2=  [<*>ML2sin(ut  )cos( <*>t  )  ] 


the  perturbation  term  becomes 


If  we  let 


d  =  M"' 


.  5<»>ML  sin(  ut  )cos((«Jt ) 

I  +  . 25ML2sin2 ( wt ) 

then 


( C-3  ) 
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*02}  [l  +  .25ML2/I02sin2(«t ) 


This  can  be  written  as 


.  P-] 

UnJ 


1  +  X) 


,  .  2.  .  .  ,  ML 

where  x  =  r  sin  ( ut  )  and  r  =  tv 
o  o  4 1 , 


(C-4  ) 


By  substituting  Eq.  (C-4)  into  (C-3)  we  get 


5wML  .  ,  .  .  „  2 

— -  sin(wt)cos(ut)( 1-  x+  x 

i02 


Expanding  further, 


(C-5  ) 


r  ;  ^ 

1 2  '  r 

— —  =  2wr  I  s  in  (  wt  )  cos  ( <*>t  )  -  sin(wt  )cos(ut  )  > 

.2J  L 

+  sinlut  )cos(ut)  x2]  ( 


(C-6  ) 


lil  . 


— —  =  2<Jr0  sin  (  wt  )  cos  ( wt  )  -  sin  (  wt  )  cos  (  wt  )  rQs  in  (  wt  ) 

2J  L 


+  s in( wt )cos ( wt )  rQsin  (wt) 


(  C-7  ) 


Using  similar  trigonometric  identities  to  those  used  in 
Appendix  B,  Eq.  (C-7)  becomes 


I  1 

~r  =  2wr0[Asi 


Asin(2ut)  +  Bsin(3wt)  +  Csin(4wt) 


+  Dsin(5wt)  +  Esin(6ut) 


( C-8  ) 


where 


A  =  1(1/2)  -  r0/4  +  (r0/2)' 
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B  =  (-3r0V32) 

C  =  (r>  -  r02/32) 

D  =  (-3r02/32) 

E  =  (r02/32) 


I  )  WT  2  ,  .  r„ 

02  ,  ML  ,  o 

Since  r  =  — - —  and  rQ  =  j=-  ,  then  rQ  becomes  rQ  =  - —  . 

[  I03j  103  2 

/ 

Because  rQ  has  a  maximum  value  of  .25,  then  all  terms  that 


are  of  higher  order  than  ^r0j  1 
Eq .  { C  —  8  )  then  becomes 


can  be  ignored. 


lz\  '  '  ro 

— —  2urQ  (  1 /2  )  sin(  2wt  )  =  WrQs  in  (  2i*>t  )  =  o  —  sin(2wt) 


(  C-9  ) 


Substituting  Eq.  (C-9)  into  Eq.  (C-2)  yields  at  this  point 


I  -  I  -  I  I  -  I 

- - -  ny  +  4  - - -  n  r  + 

i2  .2 


2  c  ^nWy  +  Wr  j  s i n ( 2wt ) 


(C-10  ) 


where  2c  = 


V  I2- 


Coef f icient 


Recalling  again  the  relationships  for  I  and  I  from 


Appendix  B,  the  coefficient  becomes 
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I  -  I  -  I 

3  2  1 


I03~  J02~  *01 _ 

1Q2  +  . 25ML2s in2 ( wt ) 


(C-l  1  ) 


I03~  I02  X01 


where  D  is  defined  in  Eq.  (C-4). 


By  defining  X2  = 


Eq.  { C  —  1 1  )  becomes 


^03  ^02  ^01 


1 ~  r2~  r 1 


1 3_  I2“  IJ  2 

— - - - 1  =  X2  (1  -  x  +  x  ) 

1  2 


(C-l 2  ) 


=  X2  -  X2  rQsin  (wt)  +  X2  rQ  sin  (ut) 
Using  trigonometric  relationships,  Eq.  (C-12)  becomes 


I  -  I  -  IJ  |  r  3 r 

-J—r —  =  >2  1_  r  +  ~r 
2 


(C-13  ) 


'  r0  r02  ]  f  r  2  ] 

+  X2  —  +  —  cos(2wt)  +  X2  g—  cos(4wt) 


Recalling  that  rp  =  —  =  2c,  Eq.  (C-13)  can  be  written  as 

1  2 


I  -  I  -  I, 

3  2  1 


+  X2 


<z)  * 12  (c  ~  i  £2) 

^  e2j  cos  (  4wt  ) 


cos ( 2ut ) 


(C-14  ) 


By  ignoring  terms  that  are  of  higher  order  than  c  ,  Eq. 
(C-14)  can  be  simplified  to 


I  -  I  -  I, 

3  2  1 


=  X2 


+  X2  e  cos ( 2wt ) 


(  C  —  1 5  ) 
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V  xi 


Coe  f  f ic ient 


This  term  can  be  written  as 


V  ul 

I„  +  . 25ML2sin2  (  wt )  -  I 

03  01 

U  J 

I  +  . 25ML2sin2 ( wt ) 

By  defining  «2 


•ik-T.") 

I  »-r,  1 


1  -  X  + 


*2) 


(C-16  ) 


L  r2 


■J 


,  En  .  (C-16)  becomes 


'  V  ^  .  («2  ♦*)(!-*♦  x2) 


(C-l 7  ) 


=  «2  +  (  1 -0(2  )  rQs  in2  ( wt  )  +  («2-l)  r^sin^fut)  +  ro3sin6(wt) 
Using  trigonometric  relationships,  Eq .  ( C  —  1 7  )  become- 

F) 


Ur  U 


„  „  3  „  2  3  2 

a2  +  e  -  «2c  +  ^*2c  -  -^c  + 


a2 e  -  e  -  2a2e2+  2c2-  ^-|c3 1  cos  (  2«t  ) 

4  ) 


+  £  |«2e2-  -|c2+  |c3jcos(4wt) 

+  ■|-E3j  cos (  6wt  ) 

Ignoring  terms  of  higher  order  than  e1  yields 

|  o(2  +  €(l-a2)|  +  I  e(a2-l  )cos(  2«t ) 


( C-18  ) 


V 


=  (o 


( C-l 9 ) 


81 


f  I  -  I  -  I  1 

3  2  1 

and 

r  i  -  i  i 

3  2 

I 

I 

l  1  J 

1  1  > 

Coef  f ic ients 


The  expansion  of  these  coefficients  yields 


A 

f  i  -  i  -  i  i 

f  .  \ 

1  -  r  -  r 

03  02  01 

2  1 

“ 

I 

r 

y 

l  01  J 

l  1  J 

’  3  2 


f  I  -  I  1 

f  l-r  1 

03  02 

2 

“ 

1 

r 

j 

l  01  J 

L  1  J 

(C-20  ) 


(C-21  ) 


As  a  result,  Eq .  (c-1)  remains  essentially  the  same  and  does 

not  contribute  to  the  excitation  of  the  system. 


\  + 


f  I  -  I  -  I  'I 

03  02  01 

nr  + 

I03-  !02  ] 

:o.  J 

:o,  J 

2 

n  y 


( C-22  ) 


Summary 

Substituting  Eqs .  ( C  —  1 5  )  and  (C-19)  into  Eq.  ( C  —  1 0  ) 


yields 


r  -  X2  1-  e 


ny  +  4  |a2  +  e(l-oc2)jn2r  + 

+  2  c  ^nwy  +  wr  jsin(2wt)  =  0 
+  2e^-  -|-X2ny  +  2  (a2-l  )n2rj  cos(  2wt ) 


( C-23  ) 
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Appendix  D 


Appendix  D  contains  the  Fortran  computer  program  used 
to  numerically  integrateEqs .  (17a),  (17b),  (17c),  (a-7a), 

(A- 7b),  (A-7c),  and  (A-7d). 
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program  non linear 


common  /ham/  t , x ( 7 , 4 ) , f ( 7 , 4 ) , er rest ( 7 ) , n , h , i 0 1 , i 02 , i 03 , m , 1 1 , 

+  w , e , sma , u , n4 

double  precision  t , x , f , errest , h , iOl , i02 , i03 , m , 11 , w , e , sma , u , n4 , 
+  z,phi,theta,psi, phi dot , thetadot ,psidot, 

+  a , check , conversion 

dimension  z(7),a(3,3) 

open! unit=8,file= ’pitch’ ,status=’new’ ) 
open(unit=9, file=’roll’ ,status=’ new ’ ) 
open(unit  =  10,file=’yaw’  ,status=’new’  ) 
open(unit  =  ll,file=’inputl’  ,status=’old’  ) 
open! unit  =  12 , file=’ input 2  *  ,status=’old’  ) 

t 1 ines=0 . 


c 

c  read  in  number  of  equations,  plotted  points, 
c  steps  between  printing,  and  timestep 

c 

•  read  (11,5)  n , npr int , nstep , t , h 

5  f ormat ( 3 1 6 , 2 f 6 . 2 ) 

write)*,*)  n , npr int , nstep , t , h 

c 

c  read  in  static  inertias,  and  m , 1 , w , e , sma , u , n4 

c 

•  read (12,*)  i01,i02,i03 

write!*,*)  i01,i02,i03 
read (12,*)  m,ll,w 
write!*,*)  m,ll,w 
read(12,*)  e,sma,u 
write!*,*)  e,sma,u 

•  n4=sqrt ( u/sma**3 ) 
c 

c  read  in  initial  conditions 
read(12,*)  theta , phi , psi 
write!*,*)  theta , phi , psi 
read! 12 ,* )  thetadot,phidot,psidot 

•  write!*,*)  thetadot , phidot , psidot 
c 

c  convert  euler  angles  to  radian  measure 
conversion  =  0.01745329252 
theta  =  theta*conversion 
phi  =  phi*conversion 

•  psi  =  psi*conversion 
c 

c  compute  initial  direction  cosine  matrix 

a(l,l)  =  cos ( phi ) *cos ( theta ) -sin ( phi ) *s in( ps i ) *s in( theta ) 
a(l,2)  =  cos ( phi ) *sin( theta ) +sin ( phi ) *s in! ps i ) *cos ( theta ) 
a(l,3)  =  -sin( phi )*cos(psi ) 

•  a(2,l)  =  -cos ( psi ) *s in( theta ) 

a!2,2)  =  coslpsi )*cos( theta  I 
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a( 2 , 3  )  =  sin( psi  ) 

a(3,l)  =  s in ( phi ) *cos ( theta ) +cos ( phi ) *s in ( ps i ) *sin ( theta ) 
a(3,2)  =  s in ( phi ) *s in ( theta ) -cos ( phi ) *s in ( ps i ) *cos ( theta ) 
a(3,3)  =  cos ( phi ) *cos ( ps i  ) 

compute  initial  angular  velocity  (  z ( 5  )  ,  z ( 6  )  ,  z ( 7  )  are  wl,w2,w3) 
z(5)  =  a ( 1 , 3 ) *n4 - thetadot*cos ( ps i ) *s in ( phi ) +ps idot*cos ( phi  ) 
z(6)  =  a ( 2 , 3 ) *n4  +  thetadot*s in ( psi  )  +phidot 

z(7)  =  a ( 3 , 3 ) *n4+ thetadot ♦cos ( psi ) *cos ( phi ) +psidot *s in ( phi ) 

compute  initial  quaternion 

z  (  4 )  =  0.5*sqrt( 1 .+a( 1 , 1 )  +  a( 2 ,2 )+a( 3,3  )  ) 

2(1)  =  ( 0 . 25/z( 4  )  )*( a( 2 , 3 )-a( 3 , 2  )  ) 
z(2)  =  (0.25/z(4) )*(a(3,l)-a(l,3)  ) 
z  (  3 )  =  ( 0 . 2  5/z ( 4  )  )*(a(l,2)-a(2,l  )  ) 

check  sum  of  quaternion  components  squared  being  equal  to  one 
check  =  z( 1  )**2  +  z( 2 )**2  +  z( 3  )**2+z( 4  )**2 

rename  variables  for  haming-su i tabl e  iteration 
x  (  1  ,  1  )  =  z  (  1  ) 
x(  2  ,  1  )  =  z(  2  ) 
x  (  3  ,  1  )  =  z  (  3  ) 
x  (  4  ,  1  )  =  z  (  4  ) 
x  (  5  ,  1  )  =  z  (  5  ) 
x  (  6  ,  1  )  =  z  (  6  ) 
x  (  7  ,  1  !  =  z  (  7  ) 

initialize  haming 
nxt  =  0 

call  haming(nxt) 
if(nxt  . ne .  0)  go  to  50 
write  ( *  ,  1  ) 

format(2x,’  haming  did  not  initialize’) 
stop 

continue 

integrate  ode.... two  nested  loops 
do  200  ipr  =  l,nprint 
do  100  istp  =  l,nstep 

each  call  to  haming  advances  one  step... 
call  haming(nxt) 


100  continue 


c  after  nstp  integration  steps,  print  current  values 

z ( 1 )  =  x ( 1 , nxt ) 
z(2)  =  x ( 2 , nxt ) 
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c 

c 

c 


c 


c 


2(3)  =  x ( 3 , nx  t ) 
z(4)  =  x ( 4 , nxt ) 

2(5)  =  x ( 5 , nxt ) 

2(6)  =  x ( 6 , nxt  ) 
z(7)  =  x ( 7 , nxt ) 

compute  direction  cosine  matrix  a,  euler  angles,  check  quaternion 

also  store  data  for  plotter 

all, 3)  =  2*( z( 1  )*z( 3  )-z( 2  )*z( 4  )  ) 
a (  2  ,  1  )  =  2*( z( 1 )*z( 2 )-z( 3  )*z( 4  )  ) 
a ( 2 , 2  )  =  -zl 1  )**2  +  z( 2 )**2-z( 3 )**2  +  z( 4  )**2 
a ( 2 , 3  )  =  2*(  z(  2 )*z( 3 )  +  z( 1  )*z( 4  )  ) 
a ( 3 , 3  )  =  -z( 1 )**2-z( 2 )**2+z( 3  )**2  +  z( 4  )**2 

theta  =  atan< -a( 2 , 1  ) /a ( 2 , 2  )  ) 

psi  =  asin( a( 2 , 3  )  ) 

phi  =  atan ( -a( 1 , 3  ) /a ( 3 , 3  )  ) 


c 


c 


c 


13 


c 

200 

c 


c 


thet  adot 


+ 

psidot  = 
phidot  = 

+ 

+ 


=  ((z(7)-a(3,3)*n4)*cos(phi)-(z(5)-a(l,3)*n4)* 
sin ( phi ) )/cos( psi  ) 

( ( z( 5 )-a( 1 ,3 )*n4 ) *cos ( phi )  +  ( z( 7 )-a( 3 , 3 )  *n4 ) *s in( phi  )  ) 
( z( 6 ) -a ( 2 , 3 ) *n4 )- 

(  ( z( 7  )-a( 3 , 3  )  *n4  )*cos( phi ) - ( z ( 5 ) -a ( 1 , 3 )*n4 )* 
sin(phi) )*tan(theta) 


check  =  z( 1 )**2  +  z( 2 )**2  +  z( 3  )**2  +  z( 4  )**2 


tlines=tlines+l . 
torbits=t/5250. 
write(8,13)  to rb its, theta 
f  ormat ( Ix,fl2.6,lx,fl2.6) 
write(9,13)  to rb its, phi 
write(10,13)  torbits,psi 
if(tlines  .ge.  14000.)  then 
close ( unit=8 ) 
close! unit  =  9  ) 
close) unit=10 ) 

open(unit=8,file=’ pitch’ ,status=’new’ ) 
open(unit=9, file=’roll’ ,status=’ new ’ ) 
open( unit=10,file=’ yaw ’ ,status=’new’ ) 
tl ines=0 . 
end  i  f 


continue 


stop 

end 


c 


subroutine  haming(nxt) 
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c 

c 

c 

c 

c. 

c 

c 


haming  is  an  ordinary  differential  equations  integrator 
it  is  a  fourth  order  predictor-corrector  algorithm 
which  means  that  it  carries  along  the  last  four 
values  of  the  state  vector,  and  extrapolates  these 
values  to  obtain  the  next  value  (the  prediction  part) 
and  then  corrects  the  extrapolated  value  to  find  a 
new  value  for  the  state  vector. 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

10 


20 


the  value  nxt  in  the  call  specifies  which  of  the  4  values 
of  the  state  vector  is  the  "next"  one. 

nxt  is  updated  by  haming  automatically,  and  is  zero  on 
the  first  call 

the  user  supplies  an  external  routine  rhs(nxt)  which 
evaluates  the  equations  of  motion 

common  /ham/  x , y ( 7 , 4 ) , f ( 7 , 4 ) , errest ( 7 ) , n , h , i 0 1 , i02 , i 03 , m , 1 1 , 
+  v , e , sma , u , n4 

double  precision  x, y , f , errest, h,hh,xo,tol,i01,i02,i03,m, 11, 

+  w , e , sma , u , n4 


all  of  the  good  stuff  is  in  this  common  block, 
x  is  the  independent  variable  (  time  ) 
y(7,4)  is  the  state  vector-  4  copies  of  it,  with  nxt 
pointing  at  the  next  one 

f ( 7 , 4  )  are  the  equations  of  motion,  again  four  copies 
a  call  to  rhs(nxt)  updates  an  entry  in  f 
errest  is  an  estimate  of  the  truncation  error  -  normally  not 
used 

n  is  the  number  of  equations  being  integrated 
h  is  the  time  step 

tol  =  1.0d-08 

switch  on  starting  algorithm  or  normal  propagation 
if(nxt)  190,10,200 

this  is  hamings  starting  algorithm. ...  a  predictor  -  corrector 
needs  4  values  of  the  state  vector,  and  you  only  have  one-  the 
initial  conditions. 

haming  uses  a  Picard  iteration  (slow  and  painfull)  to  get  the 
other  three. 

if  it  fails,  nxt  will  still  be  zero  upon  exit,  otherwise 
nxt  will  be  1,  and  you  are  all  set  to  go 

xo  =  x 

hh  =  h/2 . 0d+00 
call  rhs ( 1  ) 
do  40  1  =  2,4 
x  =  x  +  hh 
do  20  i  =  1 ,  n 

y ( i , 1 )  =  y ( i , 1-1  )  +  hh*f(i,l-l  ) 
call  rhs{ 1  ) 
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30 

40 


50 


70 


90 


110 

120 


130 

140 

150 


160 


190 


c 

c 


c 

200 


c 

210 

220 

230 


c 


x  =  x  +  hh 
do  30  i  =  1  ,  n 

y( i ,1 )  =  y ( i , 1-1  )  +  h*f ( i  ,1  ) 
call  rhs ! 1 ) 
j  sw  =  -10 
isw  =  1 

do  120  i  =  1  ,  n 

hh  =  y ( i  ,  1  )  +  h*(  9 . 0d  +  00* f ( i  ,  1  )  +  1 9 . 0d  +  00* f ( i , 2  ) 

1  -  5 . 0d  +  00*  f ( i  ,  3  )  +  f (  i  ,  4  )  )  /  24.0d  +  00 

if(  dabs(  hh  -  y(i,2))  .It.  tol  )  go  to  70 
isw  =  0 
y ( i , 2 )  =  hh 

hh  =  y  (  i  ,  1  )  +  h*{  f ( i , 1 )  +  4 . 0d+00* f ( i , 2 )  +  f ( i , 3  )  ) /3 . 0d  +  00 
if(  dabs(  hh-y(i,3))  .It.  tol  )  go  to  90 
isw  =  0 
y ( i , 3  )  =  hh 

hh  =  y ( i , 1  )  +  h* (  3 . 0d  +  00*  f ( i , 1  )  +  9 . 0d  +  00* f ( i , 2  ) 

1  +  9 . 0d  +  00*  f ( i , 3  ! 

1  +  3 . 0d  +  00*  f  {  i  ,  4  )  )  /  8 . 0d+00 

if(  dabs { hh-y ( i , 4 ) )  .It.  tol  )  go  to  110 
isw  =  0 
y ( i , 4 )  =  hh 
continue 
x  =  xo 

do  130  1  =  2,4 

x  =  x  +  h 

call  rhs (1) 

if(isw)  140,140,150 

j  sw  =  jsw  +  1 

if (jsw)  50,280,280 

x  =  xo 

isw  =  1 

jsw  =  1 

do  160  i  =  1 , n 
errest(i)  =  0.0 
nxt  =  1 
go  to  280 
jsw  =  2 

nxt  =  iabs ( nxt 1 

this  is  hamings  normal  propagation  loop  - 
x  =  x  +  h 

npl  =  mod(nxt,4)  +  1 
go  to  ( 210 , 230 ) , isw 
permute  the  index  nxt  modulo  4 
go  to  (270,270,270,220) , nxt 
isw  =  2 

nm2  =  modi npl, 4)  +  1 
nml  =  mod(nra2,4)  +  1 
npo  =  mod(nml,4)  +  1 
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c 

c 


240 


c 

c 

c 

c 


250 

260 

270 

280 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 


this  is  the  predictor  part 
do  240  i  =  1 , n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00*h*(  2 . 0d+00* f ( i , npo )  -  f(i,nml) 
1  +  2 . 0d+00*  f ( i , nm2  )  )  /  3.0d  +  00 

y(i,npl)  =  f(i,nm2  )  -  0 . 925619835d+00  *  errestl i  ) 

now  the  corrector  -  fix  up  the  extrapolated  state 
based  on  the  better  value  of  the  equations  of  motion 

call  rhs { npl ) 
do  250  i  =  1  ,  n 

v(i,npl)  =  (  9 . 0d  +  00*y (  i , npo  )  -  y(i,nm2) 

1  +  3.0d+00*h*(  f(i,npl) 

1  +  2 . 0d  +  00* f ( i , npo  )  -  f(i,nml)  )  )  /  8.0d  +  00 

errest(i)  =  f(i,nm2)  -  y(i,npl) 

y(i,npl)  =  y(i,npl)  +  0 . 074 3801 653d  +  00  *  errestl i  ) 

go  to  ( 260 , 2 70  )  , j sw 

call  rhs ( npl  ) 

nxt  =  npl 

return 

end 

subroutine  rhs(nxt) 

rhs  is  the  right  -  hand  -  side  subroutine,  customized 
to  the  particular  set  of  odes  being  integrated 

common  /ham/  t , x ( 7 , 4 ) , f ( 7 , 4 ) , errest ( 7 ) , n , h , iOl , i02 , i03 , m , 1 1 , 

+  w , e , sma , u , n4 

double  precision  t , x , f , errest , h , iO 1 , i02 , iO 3 , m , 1 1 , w , e , sma , u , 

+  n4,z,a,rl,r,nl,n2,n3,il,i2,i3,b 

dimension  z( 7 ) , r( 3 ) ,a( 3, 3  ) 

variable  description 

calculate  current  values  of  inertia 


il  = 

iOl 

i2  = 

i02  +  m*  ( ( .5*ll*sin(w*t) )**2 ) 

i  3  = 

i03  +  m*  ( ( .5*ll*sin(w*t) )**2  ) 

convert 

variable  names 

z  (  1  ) 

x ( 1 , nxt ) 

z(  2  ) 

= 

x ( 2 , nxt ) 

z  (  3  ) 

= 

x ( 3 , nxt ) 

z(  4  ) 

= 

x ( 4 , nxt ) 

z  (  5  ) 

= 

x ( 5 , nxt ) 

z  {  6  ) 

= 

x ( 6 , nxt ) 

z(  7  ) 

x { 7 , nxt ) 

calculate 

direction  cosine  matrix 

a ( 1 , 1 )  =  z( 1  )**2  -  z  (  2  ) **2  -  z ( 3  ) **2  +  z(4)**2 


89 


c 

c 

c 

c 


c 

c 


c 

c 

c 


c 


c 


c 


c 


c 


c 


a( 1 ,2  ) 
a  (  1  ,  3  ) 
a  (  2 , 1  ) 
a  (  2 , 2  ) 
a  (  2 , 3  ) 

a( 3,1  ) 
a( 3,2  ) 
a  (  3 , 3  ) 


2*( z( 1 )*z( 2  )  +  z( 3  )*z( 4  )  ) 

2* ( z ( 1  )  *z ( 3  )  -  z( 2  )*z( 4  )  ) 

2*(z(l)*z(2)  -  z( 3  )*z( 4  )  ) 

-z(l)**2  +  z{ 2  )**2  -  z ( 3  )  **2  +  z ( 4  )  **  2 
2*( z( 2  )*z( 3  )  +  z( 1  )*z( 4  )  ) 

2*( z( 1  )*z( 3  )  +  z( 2  )*z( 4  )  ) 

2*( z ( 2 )*z( 3  )  -  z( 1 )*z( 4  )  ) 

-z ( 1 ) **2  -  z ( 2  )  **2  +  z ( 3  )  **2  +  z ( 4  )  **2 


compute  rl  (length  of  orbit) 
rl  =  sma*(l  -  e*cos(n4*t)) 


calculate  radius  vector  components 
r ( 1  )  -  a ( 1 , 1  )*rl 
r( 2 )  =  a( 2 , 1  )*rl 
r( 3 )  =  a ( 3 , 1 )*rl 

calculate  gravity  gradient  torque  n 
b  =  ( 3*u ) / ( rl **5  ) 
nl  =  b*  (  i  3  —  i  2  )*r(  2  )*r(  3  ) 
n2  =  b*( i 1 - i 3 )*r( 1 )*r( 3  ) 
n3  =  b*(  12- i 1 )*r( 1 )*r(  2  ) 


write  equations  of  motion 

f(l,nxt)  =  0 . 5* ( x ( 2 , nxt ) * ( x( 7 , nxt ) ~a( 3 , 3  )  *n4  ) -x( 3 , nxt  )  * 

+  ( x( 6 , nxt ) -a( 2 , 3 ) *n4 ) +x ( 4 , nxt  )* ( x( 5 , nxt  )  -a ( 1 , 3  )  *n4  )  ) 

f(2,nxt)  =  0 . 5*  (  -x ( 1 , nxt ) * ( x ( 7 , nxt ) -a( 3 , 3 ) *n4 ) +x ( 3 , nxt ) * 
+  (x(5,nxt)-a(l,3)*n4)+x(4,nxt)*(x(6,nxt)-a(2,3)*n4)) 

f(3,nxt)  =  0 . 5* ( x ( 1 , nxt ) * ( x ( 6 , nxt ) -a ( 2 , 3 ) *n4 ) -x ( 2 , nxt )  * 

+  (x(5,nxt)-a(l,3)*n4)+x(4,nxt)*(x(7,nxt)-a(3,3)*n4)) 

f(4,nxt)  =  0 . 5* ( -x ( 1 , nxt  )  * ( x( 5 , nxt ) -a( 1 , 3 ) *n4 ) -x ( 2 , nxt ) * 
+  (x(6,nxt)-a(2,3)*n4)-x(3,nxt)*(x(7,nxt)-a(3,3)*n4)) 


f(5,nxt)  =  ( nl + ( i2- i 3 ) *x ( 6 , nxt ) *x ( 7 , nxt ) ) / i 1 

f(6,nxt)  =  ( n2+ ( i 3- i 1 ) *x ( 7 , nxt ) *x ( 5 , nxt ) - ( 0 . 5*w*m*ll**2* 
+  sin(w*t)*cos(w*t))*x(6,nxt))/i2 


f(7,nxt)  =  (  n3  +  (  i  1  -  i  2  )  *x  (  5  ,  nxt )  *x  (  6  ,  nxt  )  -  (  0 . 5*w*m*l  1**2!*: 
+  sin(w*t)*cos(w*t))*x(7,nxt))/i3 


return 

end 
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